The Northwest Atlantic Ocean one of the fastest warming locations on Earth. Research on the impacts of this rapid warming has primarily focused on high-profile and/or upper trophic level species. The temperature size rule (TSR) suggests that ectotherms living at elevated temperatures will reach smaller body-sizes at comparable development stages. However, it is unclear whether that relationship can be mitigated against at an individual level through adaptive behaviors in an open ocean environment. Here we’ve investigated ecosystem wide impacts using size spectrum analysis to track changes in the community size structure. In cases where community responses fail at adapting to elevated temperatures, we anticipated a steepening of the body-size spectrum slope relating to a reduction in larger sized individuals and an increased prevalence of smaller sized individuals. Using data from fisheries independent surveys we estimated size spectrum indices for four regions along the US NE continental shelf. Correlation/regression analyses were then performed to then assess the degree to which these changes were in alignment to hypothesized bottom-up and top-down disturbances. At the regional scale, we found that community size structure changes (spectra slope) were the largest in the Northern regions, in the Gulf of Maine and Georges Bank. These areas are home to the coldest temperatures and the largest proportions of groundfish species in the community. Spectrum slope declines were most pronounced in the 80’s and 90’s, before to the more-recent elevated temperatures. This suggests that external forces besides temperature drove the initial declines of larger-sized individuals within the communities, before elevated temperatures began to influence the ecosystem. Correlation analyses reveal that while fisheries landings are strongly correlated with these declines, bottom-up factors of zooplankton community metrics, Gulf Stream Index, and SST anomalies are also important. While the primary pressure of fisheries exploitation has declined dramatically over time, the recovery of larger-sized individuals has not been seen. That kind of recovery will likely depend on the elevated temperatures seen over the last decade.
Introduction
Temperature & Ecology
Temperature plays a critical role on biological life impacting many of the chemical reactions that underpin basic physiological function. Temperature has direct and indirect impacts on critical biological functions including the acquisition of biomass through feeding, the rates of biomass loss through metabolism, and the rates at which individuals mature and develop. Because of these relationships, most species have evolved thermal preferences around which these chemical reactions are optimized. Species that are unable to maintain their thermal preferences internally must be able to follow their thermal preference in the environment through locomotion or adapt to less-favorable conditions through changes in behavior or risk metabolic costs in failing to do so. In an era of anthropogenic climate change, there is an expectation that many species will be displaced from historic habitats in their efforts to follow their thermal preferences. Recent research in marine environments has shown evidence of this as species are now shifting to higher latitudes and to deeper depths in the pursuit of more favorable conditions ( Kleisner et al. (2017)Pinsky et al. (2013) ). Others have suggested that temperature related impacts may not be seen through geographic distribution change, but through physiological changes, changes in seasonal phenology, or in dashed hopes of species recovery ( Daan et al. (2005)Miller, O’Brien, and Fratantoni (2018)A. Pershing et al. (2018)University of South Carolina et al. (2021) ).
Need to connect temperature to size here
The potential for elevated temperatures to impact the size structure of an ecosystem has implications for the ecosystem resilience in the face of climate change, as well as the blue economies & natural resource systems that rely upon their good health.
Size Spectrum Theory
Size is a defining characteristic of all species and mediates many ecological interactions ( Brown, West, and Enquist (2000) ). Body size impacts the mobility of an organism, its ability to evade predation, its ability to find success preying on other organisms, its efforts to locate and follow essential habitats, and the metabolic costs associated with each of these behaviors. Body size also mediates how vulnerable individuals are to features their immediate environment such as temperature and/or resource availability.
Size structured environments are a fundamental organizational pattern emerging from these relationships that has been heavily researched. In the context of a strongly size-structured ecosystem, growth and maturity changes alter fitness and ultimately determine whether a species is successful in the given environment. Ecological theory is rich with models relating how energy transfers from smaller prey species to larger predatory trophic levels, the allocation of energy for growth, and the trade offs of allocating energy towards reproduction at the expense of growth. A persistent pattern in ecological systems globally is the relationship between body size and abundance.
^ Add something about GSDR to just state what the relationship is
Size spectrum and individual size distribution (ISD) models have gained increasing attention as an entry point to assessing ecosystem health and to detect systemwide disturbance ( Shin et al. (2005) ). These models avoid the need to explicitly articulate each predator-prey interactions and can be estimated from the commonly collected measures of abundance and size. A “size spectrum” describes the distribution of biomass or abundance as a function of individuals’ mass or size on a log–log scale Guiet, Poggiale, and Maury (2016) . Size spectra are described by two terms, the size spectrum slope & intercept. These two terms convey the baseline productivity, and how energy flows through an ecosystem in the form of biomass. The spectrum intercept value captures the richness or the productivity at the base of the community and is strongly connected to the prevailing environmental conditions Boudreau and Dickie (1992). So much so, that spectrum shape is sometimes defined by its eutrophic or oligotrophic environmental conditions Rossberg (2012).
Size spectra condense the complexities of predator prey networks and their interactions into a handful of size-based indices (SBI). In doing so a community size spectrum captures the emergent properties of a system, while needing only size and abundance information that is commonly available. For this reason it has become increasingly relied upon as an indicator of ecosystem health in the push for ecosystem based fisheries management. Changes in slope have been associated with fishing exploitation, primarily through the targeted removal of larger individuals ( Bianchi et al. (2000); Shin et al. (2005) ). Numerical experiments have also tried to link changes in slope to environmental disturbances Guiet, Poggiale, and Maury (2016) . Biomass spectrum present predictable intercepts between ecosystems of similar productivity levels, but also of distinct temperatures Guiet, Poggiale, and Maury (2016).
Direct quote from Guiet, Poggiale, and Maury (2016) , but nails the connection back to temp expectations:
Because it controls chemical reactions, temperature controls metabolic rates which underpin maintenance, growth or reproduction (Clarke and Johnston, 1999; Kooijman, 2010) as well as the functional responses to food density (Rall et al., 2012). Guiet et al. (2016)… In addition to the impact of temperature on communities’ intercepts (heights), the impact of temperature on the speed of the energy flow within communities may affect other properties, such as their resilience to perturbations or the intensity of trophic cascades (Andersen and Pedersen, 2009).
Species Trends in the Northeast Atlantic Shelf
The continental shelf groundfish community in Northwest Atlantic has changed dramatically over the last century. Stocks that supported international fishing effort collapsed, and recovery efforts fell short of their objectives. Research on Georges Bank estimated that biomass more than halved in the 1960’s (pre-dating federal monitoring efforts), and noted a species replacement of commercial groundfish target species by skate and dogfish Fogarty and Murawski (1998). Industrial fishing is inherently size-selective, with larger individuals selectively removed from the population. This has an immediate impact on the community size-distribution with additional impacts on the future population as well. Larger individuals have a greater impact on population recovery, capable of holding more (and often of higher quality) eggs. Size-based harvest in fisheries has been shown to create selective pressures that promote characteristics of early maturation at smaller sizes.
Temperature of the Gulf of Maine & NE Shelf
In addition to the ecological disturbance of industrial fisheries, the Northwest Atlantic is also one of the fastest warming locations in the global oceans. Sea surface temperatures in the Gulf of Maine since 1982 have been warming at rates faster than 96% of the world’s oceans, with similar warming rates along the northwest Atlantic shelf A. Pershing et al. (2018). This persistent elevated temperature regime of the area is a result of several forces, a combination of shifting ocean currents and the unique bathymetry of the region. A Northward shift in the Gulf Stream directly increased the regional temperatures through increased transport of warm Gulf Stream water into areas like the Gulf of Maine. The Northward Gulf Stream shift is associated with a higher frequency of warm core rings, and the obstruction of cold-water Scotian Shelf current flow that would otherwise counter the influence of the Gulf Stream on the region’s temperatures ( Gangopadhyay et al. (2019)University of South Carolina et al. (2021) ). The combination of these oceanographic changes has led to a warmer continental shelf habitat.
The rapid warming in the northwest Atlantic is a major factor in the redistribution of marine species along the US east coast. Species have responded by adjusting the timing and locations of their seasonal migrations and shifting their geographic ranges ( Nye et al. (2009)Staudinger et al. (2019) ). There is evidence that warming has hampered fisheries recoveries as well. Adding a metabolic tax to physiological pathways like growth and metabolism. Species like Black Sea Bass, Atlantic shortfin squid, and Blue crab have been high-profile examples of species expanding their ranges to follow their thermal preferences. While species like the American lobster have shown declines at their southern range near Long Island Sound, with much doubt whether they will recover under the present temperature trends. The recent regime shift in the physical oceanography has also shown to be a catalyst for biological shifts as well ( University of South Carolina et al. (2021) ; Perretti et al. (2017) ).
While these examples show that species can respond to changes in the physical environment around them through movement & behavior, research elsewhere suggests that physiological responses integrated across species will manifest as changes in community size structure.
Purpose
With the understanding that populations depend on the health of their ecosystems, there is a need to have community-wide metrics to effectively understand and manage marine resources. Size based indices are metrics that can be estimated from the information that has historically been available from long-term survey efforts. These indices have been shown to be sensitive to the impacts of fishing, but should also capture environmentally driven stress as well. We estimated size spectrum relationships as SBI’s for the groundfish populations for each sub-region of the Northeast US continental shelf. In the case of the NW Atlantic sustained increases in temperature should have a physiological impact on the community size structure.
This leads to our second hypothesis:
H2. Warming alters the community through the direct influence of temperature on metabolism, growth, and population productivity.
Methods
Groundfish Data
Fishery Independent data on was collected as part of the NEFSC’s northeast trawl survey. This survey is conducted each year in the spring and in the fall, with sample locations determined following a stratified-random survey design with effort allocated in proportion to stratum area. Trawls are performed for a fixed duration at each station, reporting abundance at length and total biomass for all species caught. Trawl survey analyses incorporated data from both the Spring and Fall survey seasons, and for all years from 1970 to 2019. Correction factors were applied to total species abundance and aggregate species biomass to account for all changes in vessels, gear, and doors when appropriate as part of the survey program. However, abundance and biomass at length is not corrected for as part of the standard survey data protocol and needed to be estimated. To account for this, abundance at length for each species were adjusted to match the correction factors applied to total species abundance at each station, with allocations following the distributions of length caught at that station. Such that for each species: the sum of the resulting adjusted abundance numbers across each length is equal to the total abundance that was corrected for changes in vessels, gear, etc. To account for differences in sampling effort among survey strata, all corrected abundance-at-length data was then area-stratified.
Data from the survey was grouped by strata to form geographically meaningful sub-regions: Gulf of Maine, Georges Bank, Southern New England, Mid-Atlantic Bight. For each region, we developed several time series indicators:
Community Composition metrics (abundance and biomass by functional group, with body-size contributions)
Mean size of the aggregate community and key functional groups
Slope and intercept of the size spectrum
Community Composition
Functional groups were assigned to each species based on life history and geography using the definitions of Hare et al. (2010). Functional groups included were coastal, diadromous, elasmobranch, groundfish, pelagic, and reef species. Stratified annual abundance and biomass totals were calculated for each functional group and each region with labels for increasing body-size (biomass kg) groups.
Body Size Trends
Annual abundance weighted body length and body weight within each region and for each functional group were also estimated using the numbers at length and the estimated biomass at length information. Data for body size trends were not truncated and include all sizes caught.
Size Spectrum Analyses
Community size spectra were estimated using abundance-at-length data from 68 species. These species were selected based on the availability of published weight-at-length relationships (Wigley et al. 2003) and represented 98.98% of the total biomass caught in the survey. Published length-weight relationships were used to convert abundance at length data into their corresponding biomass at length (kg). These values were then used to get totals for stratified weight-at-length, in complement to the corrected abundances-at-length data which had been area stratified. These area-stratified biomass at length totals were then used for fitting each regional biomass size spectra and exponents of individual size distributions. Data for these analyses was truncated by body masses of a minimum 1g and a maximum 10kg to account for poor gear selectivity at the smallest and largest size ranges and to ensure that all size bins were sampled in all years for the body mass spectrum fitting.
Body Mass Spectra
Normalized biomass spectra were estimated following the LBNbiom methodology as described in A. M. Edwards et al. (2017). When fitting the normalized biomass size spectra, stratified biomass at length data was binned into equal spaced intervals on a (1) on a \(log_{10}\) scale, with bodymass totaled across all species. To normalize the spectra, the stratified abundances within each bin was then divided by the bin-width to account for the increasing bin-widths, a consequence of the log scale. Normalized bodymass spectra were fit for each year and for each region independently, and for each year across all strata, using ordinary least squares (OLS) regressions for stratified abundance (normalized) by body-size bins.
Individual Size Distributions
Individuals in the catch data are measured to the nearest cm, with smaller individuals measured to the nearest millimeter. Because individual biomass is estimated from those length measurements, there is for each increment a range of possible bodymass that spans between the cm & mm measurement increments. The relationships between length and bodymass in fishes is exponential and taxon specific, so biases resulting from using only the lower or upper end of those ranges is different for each taxon and increases for larger taxa. To account for this, we used the extended likelihood method (MLEbin) of A. Edwards et al. (2020) . This method estimates the exponent of size spectra (b) for a bounded power law relationship between individual abundances and their estimated biomass. With the MLEbin method the ranges of bodymass associated with each length measurement is accounted for to reduce bias. Exponents of the individual size distribution were estimated for each region and for every year. A minimum body mass of 1g was used for the lower bound with 10kg individual bodymass used as an upper bound for the ISD relationship. The biomass within these limits represents 97.83% of all estimated biomass used in these analyses.
Temperature Data
Global Sea surface temperature data was obtained via NOAA’s optimally interpolated SST analysis (OISSTv2), providing daily temperature values at a 0.25° latitude x 0.25° longitude resolution (Reynolds et al. 2007). A daily climatology for every 0.25° pixel in the global data set was created using average daily temperatures spanning the period of 1982-2011. Daily anomalies were then computed as the difference between observed temperatures and the daily climatological average. OISSTv2 data used in these analyses were provided by the NOAA PSL, Boulder, Colorado, USA from their website at https://psl.noaa.gov.
Temperature data was regionally averaged to match the survey regions from the age-at-length data. SST anomalies were averaged by year for each region and over the entire sampling region to produce daily time series. These time series were then processed into annual timeseries of surface temperatures and anomalies. All region-averaging was done with area-weighting of the latitude/longitude grid cells to account for differences in cell-size in of the OISSTv2 data.
Spectra Drivers
The impact of external factors on the changes in size spectra was correlated against several hypothesized driving forces related to both environmental regimes and anthropogenic disturbances. Potential environmental drivers include sea surface temperature anomalies, Gulf Stream Index (GSI), and two zooplankton community indices derived from the Gulf of Maine continuous plankton recorder (CPR) program. Hypothesized anthropogenic drivers include state and federal fisheries landings from the Greater Atlantic Regional Fisheries Office (GARFO), divided by reporting zones into aggregate regions to closely align with the survey areas we defined for the size spectra analyses.
Results
Abundance Distribution
Abundance across all body sizes remained relatively stable from the 1970’s before rising in the northern regions around 1990 beginning in the Gulf of Maine. Around this time abundances increased through the mid 2010’s. Further south in Georges Bank, abundances remained flat until the 2010’s, when overall abundance roughly tripled, only to fall back to previous amounts by the end of the century. Southern New England saw a less dramatic rise and fall that began just before 2010, again falling back to earlier levels by the end of the century. The Mid Atlantic Bight had relatively consistent abundances throughout, with no major periods of abundance growth or decline, but with larger inter-annual variability.
Abundance gains observed in Georges Bank and Gulf of Maine were primarily from groundfish species, with additional growth from demersal species in the Gulf of Maine. Increases in abundance across all areas was mostly confined to individuals weighing less than .5kg. With some years driven in large-part by exceptional year-classes in a handful of species like haddock in Georges Bank. The observed abundance volatility in Southern New England and the Mid-Atlantic Bight conversely was largely the result of changes in abundance in pelagic species, whose abundance varied by several times the magnitude that of the other functional groups.
Biomass Distribution
Overall biomass was highest in the two northern regions, the Gulf of Maine and Georges Bank. Roughly half of the biomass sampled in these regions can be attributed to groundfish/demersal species, with the second largest contributions coming from elasmobranchs. Groundfish biomass, larger individuals >2kg in particular, declined during the 70’s and 80’s in these regions, never truly recovering. Beginning in the 2000’s there were signs that groundfish abundances were increasing as evidenced by increasing numbers of smaller individuals, however in both regions this trend appears to have reversed by the mid 2010’s. Elasmobranch biomass increased steadily throughout the survey time period across all regions, with the exception of southern New England. This area showed large 5-10 year swings in biomass, but no clear long-term trend. Larger elasmobranch were rare in all regions except for a period spanning the late 70’s through the early 90’s isolated to Georges Bank. Demersal species biomass was highest in the Gulf of Maine, dwarfing their contributions in other regions. Their biomass declined in the 70’s, was flat until the late 90’s, remaining relatively high until declining in the late 2010’s. Pelagic species biomass was low in all regions, and is unlikely to be representative of true biomass trends due to gear selectivity.
Regional Variation in Species Composition
There was a distinct difference between Northern and Southern regions in the way biomass was distributed among the different functional groups. The primary contributors to overall biomass in the southern regions (southern New England & mid-Atlantic bight) was the elasmobranch community. While the northern regions (Gulf of Maine & Georges Bank) each had similar quantities of elasmobranch biomasses, there was also a comparable contribution of groundfish and in the Gulf of Maine there was a major component of demersal species as well.
Body Size Trends
The average fish size in the Gulf of Maine (length and weight) declined the greatest of all regions over our study period. The average individual length was greatest in the 1970’s in the 35-40cm range, falling to 28-33cm over the last decade. Body-weight fell dramatically in the 1980’s, from around .75kg in the 1970’s to .25-.30kg, roughly a third of what it had been. Georges Bank body sizes also declined during the study period, but less dramatically. Both of these Northern regions had brief period in the late 2000’s where average length and weight rose, before falling again in the 2010’s. The MAB region was the only region to see a long-term increase in both length and weight during the study period. SNE saw no long-term change in length, and a minor decline in average body-weight.
Beginning in the 1970’s, there was a clear difference in the relative positions of spectra parameters among the different regions. Gulf of Maine and Georges Bank showed the least steep spectra slopes in the earlier time periods with slopes around -1 & -1.1 respectively. The relatively flat slopes in these regions both steepened over time, settling near -1.3 (GoM) and -1.5 (GB). Gulf of Maine experienced much of its decline during the 1980’s and 1990’s. There was a brief reversal in this trend during the 2000’s, but slopes continued to steepen by 2010 and remained steep through 2019. Georges Bank did not experience as rapid of a decline, but experienced a similar long-term steepening. In contrast to the northern regions, SNE and MAB had steeper slopes in the -1.2 to -1.5 territory. The long term pattern for SNE was one of increasing volatility, but not so much a decline. The spectra slope for the MAB was less volatile, but similarly maintained a relatively stable wander around -1.4. By the end of the study period all regions had slopes that were at or near a similar level.
Individual Size Distribution
Size Spectra Drivers
Driver Correlations
NOTE: Correlation matrix is computed starting at the year where there are no NA values in any drivers. Currently with SST included that begins the matrix at 1982.
Top-down and bottom up influences on both carrying capacity (intercept) and transfer efficiency (slope)
Some of the major drivers suggested here operate on both, but to varying degrees. Here are some potential mechanisms:
Literature suggests: - Intercept (a proxy for productivity and carrying capacity) is primarily determined by bottom up features like: nutrient availability, temperature
Slope (a measure of energy transfer efficiency and static biomass distribution) has been shown to be sensitive to the physical removal of species through fishing.
Temperature Mechanisms: - Temperature’s impact on growth via genetic plasticity impacts both the available biomass at the primary producer level (impacting ecosystem carrying capacity), as well as the Linf of larger species (recruitment rate). - Temperature also impacts the efficiency of energy (as biomass) being transferred between individuals via predator & prey interactions. More energy per-capita is expended in the form of increased metabolic rates and/or behavioral changes. This metabolic tax should steepen the spectrum slope by removing available energy at a system wide level. - Temperature impacts behavior via physiological impacts on metabolism and foraging rates as well as through the avoidance of temperature stresses.
Separating Complimentary Forces Impacting Growth
The data used for this analysis was collected as part of a survey program which began out of concern that fisheries were already being over-harvested. Early estimates from scientists at that time suggested that by the 1970’s total biomass of Georges Bank had already been halved, and elasmobranch species had begun to replace the over-harvested gadoid species (Fogarty and Murawski 1998). Having such a large disturbance which pre-dates our time series is suggestive that the measured steepening of size spectrum slope we observed in this area and the adjacent Gulf of Maine are potentially the tail-ends of a longer and more severe ecosystem decline. While metrics of overall fishing pressure do not align exactly with trawl survey coverage, historical records and anecdotal evidence fro that time suggest that groundfish fishing pressure in these areas are a fraction of their what their impacts were in the 1960’s and 1970’s.
Forces Preventing the Recovery of Large Individuals
This begs the question of why larger adult numbers never began to recover in these regions. Looking at abundance and biomass information from the survey there was evidence of strong recruitment among smaller individuals < 1kg, but there has since not been a recovery in fishes larger than 1kg outside of the elasmobranchs. Work by (A. J. Pershing et al. 2015) suggested this prolonged recovery period may be due to a lack of accounting for temperature change in fisheries management. At the time of that research, the regional temperatures had only begun rising, and could have been considered at that time an acute stressor. Since that time the region has experienced nearly a decade of sustained above-average temperatures. There are signs that the success seen in recruitment and survival of even the smaller size classes is declining. While temperature change has been associated with changes in growth rates and size-at-age, so too have size-selective fishing practices, making it difficult to disentangle the importance of exploitation & temperature on the overall community size structure when body size integrates these two forces (Shackell and Frank 2007).
Our ability to make statements on community size-structure change in this study is limited by the community that is effectively sampled as part of the groundfish survey, and of species with available weight-at-length information. When looking at the species compositions in each region it is clear that in SNE and MAB a larger proportion of the community biomass was held within a relatively small number of species, primarily the elasmobranchs. Whether or not this is incomplete representation of the community due to the sampling gear, or an accurate representation of the entire community is unclear. In these two regions the biomass spectra information is largely that of the elasmobranch community. Whether or not this is a sufficient fraction of the broader community, what key community members are absent, and what the ecological implications of this are remain unclear.
Looking at both the level of noise in the body-size trends and the size spectra parameters, the Gulf of Maine and Georges Bank seem to have the clearest signals coming through. It is likely not coincidental that these two regions also have large proportions of their sampled biomass represented by groundfish and demersal species, functional groups that are best sampled by a bottom trawl survey. Looking specifically at these two regions only, we see two parallel trends of declining body size (length and weight) among the community as a whole. We also simultaneously see a steepening of the size spectra slope. Both of these trends are slightly more severe in the Gulf of Maine.
Functional Group Assignments and Regional Presence/Absence
Common Name
Gulf of Maine
Georges Bank
Southern New England
Mid-Atlantic Bight
Coastal
Atlantic Croaker
X
X
X
Atlantic Thread Herring
X
X
Blueback Herring
X
X
X
X
Bluefish
X
X
X
X
Butterfish
X
X
X
X
Northern Kingfish
X
X
X
Southern Kingfish
X
Spanish Mackerel
X
Spanish Sardine
X
Spot
X
X
Striped Bass
X
X
X
X
Weakfish
X
X
X
Elasmobranch
Atlantic Angel Shark
X
Atlantic Sharpnose Shark
X
Barndoor Skate
X
X
X
X
Bullnose Ray
X
X
Chain Dogfish
X
X
X
Clearnose Skate
X
X
Cownose Ray
X
Little Skate
X
X
X
X
Rosette Skate
X
X
X
X
Roughtail Stingray
X
Sand Tiger
X
Sandbar Shark
X
X
Smooth Butterfly Ray
X
Smooth Dogfish
X
X
X
X
Smooth Skate
X
X
X
X
Spiny Butterfly Ray
X
Spiny Dogfish
X
X
X
X
Thorny Skate
X
X
X
X
Winter Skate
X
X
X
X
Groundfish
Acadian Redfish
X
X
X
X
American Plaice
X
X
X
X
Atlantic Cod
X
X
X
X
Atlantic Halibut
X
X
X
Atlantic Wolffish
X
X
X
Cusk
X
X
X
X
Fawn Cusk-Eel
X
X
X
X
Fourspot Flounder
X
X
X
X
Goosefish
X
X
X
X
Haddock
X
X
X
X
Longhorn Sculpin
X
X
X
X
Northern Searobin
X
X
X
X
Ocean Pout
X
X
X
X
Offshore Hake
X
X
X
X
Pollock
X
X
X
X
Red Hake
X
X
X
X
Sea Raven
X
X
X
X
Silver Hake
X
X
X
X
Spotted Hake
X
X
X
X
Summer Flounder
X
X
X
X
White Hake
X
X
X
X
Windowpane Flounder
X
X
X
X
Winter Flounder
X
X
X
X
Witch Flounder
X
X
X
X
Yellowtail Flounder
X
X
X
X
Pelagic
Atlantic Herring
X
X
X
X
Atlantic Mackerel
X
X
X
X
Buckler Dory
X
X
X
X
Round Herring
X
X
X
X
Reef
Atlantic Spadefish
X
Black Sea Bass
X
X
X
X
Blackbelly Rosefish
X
X
X
X
Cunner
X
X
X
X
Greater Amberjack
X
X
Scup
X
X
X
X
NA
American Shad
X
X
X
X
Atlantic Sturgeon
X
Functional group assignments adapted from ________
Largest Commercial Fisheries Landings of Northeastern US (by weight)
Harvest Region
Common Name
Avg. Annual Landings (lb.)
Total Landings (lb.)
Total Value ($)
1964 - 1995
Gulf of Maine
Herring, Atlantic
22.03M
1.89B
90.80M
Gulf of Maine
Menhaden, Atlantic
16.74M
954.35M
25.26M
Gulf of Maine
Hake, Silver
4.43M
566.87M
48.39M
Gulf of Maine
Cod, Atlantic
2.57M
545.39M
244.17M
Gulf of Maine
Pollock
1.91M
399.19M
107.17M
Georges Bank
Cod, Atlantic
7.67M
866.96M
411.09M
Georges Bank
Haddock
3.99M
450.77M
145.62M
Georges Bank
Flounder, Yellowtail
2.72M
307.27M
147.00M
Georges Bank
Hake, Silver
2.37M
246.08M
23.00M
Georges Bank
Flounder, Winter
1.99M
226.42M
156.36M
Southern New England
Other Fish, Bony
6.28M
627.79M
6.54M
Southern New England
Menhaden, Atlantic
6.47M
523.77M
19.26M
Southern New England
Flounder, Yellowtail
2.27M
516.96M
171.50M
Southern New England
Hake, Silver
1.88M
425.64M
98.62M
Southern New England
Flounder, Winter
1.13M
256.20M
97.99M
Mid-Atlantic Bight
Menhaden, Atlantic
74.13M
3.26B
154.62M
Mid-Atlantic Bight
Flounder, Summer
818.13K
115.36M
106.29M
Mid-Atlantic Bight
Mackerel, Atlantic
755.21K
79.30M
8.82M
Mid-Atlantic Bight
Scup
489.82K
57.80M
23.82M
Mid-Atlantic Bight
Weakfish/Sea Trout, Squeteague
446.27K
55.34M
23.15M
1996 - 2021
Gulf of Maine
Herring, Atlantic
12.26M
404.62M
24.36M
Gulf of Maine
Shark, Dogfish, Spiny
1.68M
68.84M
11.49M
Gulf of Maine
Cod, Atlantic
730.89K
65.05M
81.46M
Gulf of Maine
Monkfish/Angler/Goosefish
551.36K
49.62M
78.13M
Gulf of Maine
Pollock
593.95K
48.11M
38.27M
Georges Bank
Cod, Atlantic
2.24M
114.02M
133.29M
Georges Bank
Herring, Atlantic
4.24M
97.43M
6.54M
Georges Bank
Hake, Silver
1.43M
76.96M
32.48M
Georges Bank
Flounder, Winter
1.15M
60.01M
66.94M
Georges Bank
Haddock
1.06M
53.82M
70.62M
Southern New England
Mackerel, Atlantic
1.79M
191.14M
28.09M
Southern New England
Herring, Atlantic
2.23M
151.97M
10.92M
Southern New England
Hake, Silver
1.20M
142.41M
66.15M
Southern New England
Menhaden, Atlantic
2.26M
106.37M
8.48M
Southern New England
Skate, Nk
1.02M
101.88M
14.21M
Mid-Atlantic Bight
Menhaden, Atlantic
83.91M
6.63B
468.26M
Mid-Atlantic Bight
Croaker, Atlantic
1.81M
188.04M
82.10M
Mid-Atlantic Bight
Mackerel, Atlantic
1.73M
121.11M
14.99M
Mid-Atlantic Bight
Bass, Striped
1.30M
79.47M
173.06M
Mid-Atlantic Bight
Spot
835.86K
57.67M
41.54M
Landings data obtained from the Greater Atlantic Regional Fishing Office (GARFO)
References
Bianchi, G., H. Gislason, K. Graham, L. Hill, X. Jin, K. Koranteng, S. Manickchand-Heileman, et al. 2000. “Impact of Fishing on Size Composition and Diversity of Demersal Fish Communities.”ICES Journal of Marine Science 57 (3): 558–71. https://doi.org/10.1006/jmsc.2000.0727.
Boudreau, P. R., and L. M. Dickie. 1992. “Biomass Spectra of Aquatic Ecosystems in Relation to Fisheries Yield.”Canadian Journal of Fisheries and Aquatic Sciences 49 (8): 1528–38. https://doi.org/10.1139/f92-169.
Brown, James, G. B. West, and Brian Enquist. 2000. “Scaling in Biology: Patterns and Processes, Causes and Consequences.” In, 1–24.
Daan, Niels, Henrik Gislason, John G. Pope, and Jake C. Rice. 2005. “Changes in the North Sea Fish Community: Evidence of Indirect Effects of Fishing?”ICES Journal of Marine Science 62 (2): 177–88. https://doi.org/10.1016/j.icesjms.2004.08.020.
Edwards, Am, Jpw Robinson, Jl Blanchard, Jk Baum, and Mj Plank. 2020. “Accounting for the Bin Structure of Data Removes Bias When Fitting Size Spectra.”Marine Ecology Progress Series 636 (February): 19–33. https://doi.org/10.3354/meps13230.
Edwards, Andrew M., James P. W. Robinson, Michael J. Plank, Julia K. Baum, and Julia L. Blanchard. 2017. “Testing and Recommending Methods for Fitting Size Spectra to Data.”Methods in Ecology and Evolution 8 (1): 57–67. https://doi.org/10.1111/2041-210X.12641.
Fogarty, Michael J., and Steven A. Murawski. 1998. “Large-Scale Disturbance and the Structure of Marine Systems: Fishery Impacts on Georges Bank.”Ecological Applications 8 (sp1): S6–22. https://doi.org/10.1890/1051-0761(1998)8[S6:LDATSO]2.0.CO;2.
Gangopadhyay, Avijit, Glen Gawarkiewicz, E. Nishchitha S. Silva, M. Monim, and Jenifer Clark. 2019. “An Observed Regime Shift in the Formation of Warm Core Rings from the Gulf Stream.”Scientific Reports 9 (1): 12319. https://doi.org/10.1038/s41598-019-48661-9.
Guiet, Jérôme, Jean-Christophe Poggiale, and Olivier Maury. 2016. “Modelling the Community Size-Spectrum: Recent Developments and New Directions.”Ecological Modelling 337 (October): 4–14. https://doi.org/10.1016/j.ecolmodel.2016.05.015.
Hare, Jonathan A., Michael A. Alexander, Michael J. Fogarty, Erik H. Williams, and James D. Scott. 2010. “Forecasting the Dynamics of a Coastal Fishery Species Using a Coupled Climatepopulation Model.”Ecological Applications 20 (2): 452–64. https://doi.org/10.1890/08-1863.1.
Kleisner, Kristin M., Michael J. Fogarty, Sally McGee, Jonathan A. Hare, Skye Moret, Charles T. Perretti, and Vincent S. Saba. 2017. “Marine Species Distribution Shifts on the U.S. Northeast Continental Shelf Under Continued Ocean Warming.”Progress in Oceanography 153 (April): 24–36. https://doi.org/10.1016/j.pocean.2017.04.001.
Miller, Timothy J., Loretta O’Brien, and Paula S. Fratantoni. 2018. “Temporal and Environmental Variation in Growth and Maturity and Effects on Management Reference Points of Georges Bank Atlantic Cod.”Canadian Journal of Fisheries and Aquatic Sciences 75 (12): 2159–71. https://doi.org/10.1139/cjfas-2017-0124.
Nye, Ja, Js Link, Ja Hare, and Wj Overholtz. 2009. “Changing Spatial Distribution of Fish Stocks in Relation to Climate and Population Size on the Northeast United States Continental Shelf.”Marine Ecology Progress Series 393 (October): 111–29. https://doi.org/10.3354/meps08220.
Perretti, Ct, Mj Fogarty, Kd Friedland, Ja Hare, Sm Lucey, Rs McBride, Tj Miller, et al. 2017. “Regime Shifts in Fish Recruitment on the Northeast US Continental Shelf.”Marine Ecology Progress Series 574 (July): 1–11. https://doi.org/10.3354/meps12183.
Pershing, Andrew J., Michael A. Alexander, Christina M. Hernandez, Lisa A. Kerr, Arnault Le Bris, Katherine E. Mills, Janet A. Nye, et al. 2015. “Slow Adaptation in the Face of Rapid Warming Leads to Collapse of the Gulf of Maine Cod Fishery.”Science 350 (6262): 809–12. https://doi.org/10.1126/science.aac9819.
Pershing, Andrew, Katherine Mills, Alexa Dayton, Bradley Franklin, and Brian Kennedy. 2018. “Evidence for Adaptation from the 2016 Marine Heatwave in the Northwest Atlantic Ocean.”Oceanography 31 (2). https://doi.org/10.5670/oceanog.2018.213.
Pinsky, Malin L., Boris Worm, Michael J. Fogarty, Jorge L. Sarmiento, and Simon A. Levin. 2013. “Marine Taxa Track Local Climate Velocities.”Science 341 (6151): 1239–42. https://doi.org/10.1126/science.1239352.
Rossberg, Axel G. 2012. “6 - A Complete Analytic Theory for Structure and Dynamics of Populations and Communities Spanning Wide Ranges in Body Size.” In, edited by Ute Jacob and Guy Woodward, 46:427–521. Global Change in Multispecies Systems Part 1. Academic Press. https://doi.org/10.1016/B978-0-12-396992-7.00008-3.
Shackell, Nl, and Kt Frank. 2007. “Compensation in Exploited Marine Fish Communities on the Scotian Shelf, Canada.”Marine Ecology Progress Series 336 (April): 235–47. https://doi.org/10.3354/meps336235.
Shin, Yunne-Jai, Marie-Joëlle Rochet, Simon Jennings, John G. Field, and Henrik Gislason. 2005. “Using Size-Based Indicators to Evaluate the Ecosystem Effects of Fishing.”ICES Journal of Marine Science 62 (3): 384–96. https://doi.org/10.1016/j.icesjms.2005.01.004.
Staudinger, Michelle D., Katherine E. Mills, Karen Stamieszkin, Nicholas R. Record, Christine A. Hudak, Andrew Allyn, Antony Diamond, et al. 2019. “It’s about Time: A Synthesis of Changing Phenology in the Gulf of Maine Ecosystem.”Fisheries Oceanography 28 (5): 532–66. https://doi.org/10.1111/fog.12429.
University of South Carolina, Erin Meyer-Gutbrod, Charles Greene, Kimberley Davies, and David Johns. 2021. “Ocean Regime Shift Is Driving Collapse of the North Atlantic Right Whale Population.”Oceanography 34 (3): 22–31. https://doi.org/10.5670/oceanog.2021.308.
Source Code
---title: "Rising Temperatures Extend Community Size-Structure Damage Done by Decades of Size-Selective Fishing"author: name: "Adam A. Kemberling" title: "Quantitative Research Associate" url: https://github.com/adamkemberling affiliation: Gulf of Maine Research Institutedescription: | Results from a Size Spectrum Analysis of the Northeast US Groundfish Communitydate: "`r Sys.Date()`"# format: docxformat: html: self-contained: true code-fold: true code-tools: true df-print: kable toc: true toc-depth: 2editor: sourceexecute: echo: false warning: false message: false fig.height: 6 fig.width: 6 fig.align: "center" comment: ""bibliography: references.bib---```{r}#| label: load packages and functions#### Packages ####library(targets)library(rnaturalearth)library(here)library(sf)library(gmRi)library(patchwork)library(gt)library(knitr)library(tidyverse)library(ggstream)library(ggforce)library(bcp)library(scales)library(corrplot)library(readxl)library(ggforce)library(Hmisc)# Support functionssource(here("R/support/sizeSpectra_support.R"))# Set themetheme_set(theme_gmri() +theme(panel.grid.major =element_line(size =0.5, linetype =1, color ="gray"),panel.grid.major.x =element_line(size =0.5, linetype =3, color ="gray"),panel.grid.minor =element_line(size =0.5, linetype =3, color ="gray90"),axis.line =element_line(color ="black"),axis.line.y =element_line(color ="black"), strip.background =element_rect(colour =gmri_cols("teal")),legend.position ="bottom") )# Map polygonsus_poly <-ne_states("united states of america", returnclass ="sf")canada <-ne_states("canada", returnclass ="sf")# levels for faceting areasarea_levels <-c("Northeast Shelf, Full", "GoM", "GB", "SNE", "MAB")area_levels_long <-c("Northeast Shelf, Full", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")``````{r}#| label: group-size-summaries# Function to process summaries for various factor combinationsget_group_summaries <-function(data_in, ...){# Do some grouping to get totals group_totals <- data_in %>%group_by(...) %>%summarise(total_survey_catch =sum(numlen, na.rm = T),total_lw_bio =sum(sum_weight_kg, na.rm = T),total_strat_abund =sum(strat_total_abund_s, na.rm = T),total_strat_lw_bio =sum(strat_total_lwbio_s, na.rm = T), .groups ="drop") # length bins group_lengths <- data_in %>%group_by(..., length_bin) %>%summarise(lenbin_survey_catch =sum(numlen),lenbin_lw_bio =sum(sum_weight_kg),lenbin_strat_abund =sum(strat_total_abund_s),lenbin_strat_lw_bio =sum(strat_total_lwbio_s), .groups ="drop") %>%left_join(group_totals) %>%mutate(perc_total_catch = (lenbin_survey_catch - total_survey_catch) *100,perc_lw_bio = (lenbin_lw_bio - total_lw_bio) *100,perc_strat_abund = (lenbin_strat_abund - total_strat_abund) *100,perc_strat_lw_bio = (lenbin_strat_lw_bio - total_strat_lw_bio) *100)# weight bins group_weights <- data_in %>%group_by(..., weight_bin) %>%summarise(wtbin_survey_catch =sum(numlen),wtbin_lw_bio =sum(sum_weight_kg),wtbin_strat_abund =sum(strat_total_abund_s),wtbin_strat_lw_bio =sum(strat_total_lwbio_s),.groups ="drop") %>%left_join(group_totals) %>%mutate(perc_total_catch = (wtbin_survey_catch / total_survey_catch) *100,perc_lw_bio = (wtbin_lw_bio / total_lw_bio) *100,perc_strat_abund = (wtbin_strat_abund / total_strat_abund) *100,perc_strat_lw_bio = (wtbin_strat_lw_bio / total_strat_lw_bio) *100)return(list("length_bins"=drop_na(group_lengths),"weight_bins"=drop_na(group_weights)))}# Function to process summaries for various factor combinationsl10_bin_summaries <-function(data_in, ...){# Do some grouping to get totals group_totals <- data_in %>%group_by(...) %>%summarise(total_survey_catch =sum(numlen, na.rm = T),total_lw_bio =sum(sum_weight_kg, na.rm = T),total_strat_abund =sum(strat_total_abund_s, na.rm = T),total_strat_lw_bio =sum(strat_total_lwbio_s, na.rm = T), .groups ="drop") # weight bins group_weights <- data_in %>%mutate(weight_bin = bin_label) %>%group_by(..., left_lim) %>%summarise(wtbin_survey_catch =sum(numlen),wtbin_lw_bio =sum(sum_weight_kg),wtbin_strat_abund =sum(strat_total_abund_s),wtbin_strat_lw_bio =sum(strat_total_lwbio_s),.groups ="drop") %>%left_join(group_totals) %>%mutate(perc_total_catch = (wtbin_survey_catch / total_survey_catch) *100,perc_lw_bio = (wtbin_lw_bio / total_lw_bio) *100,perc_strat_abund = (wtbin_strat_abund / total_strat_abund) *100,perc_strat_lw_bio = (wtbin_strat_lw_bio / total_strat_lw_bio) *100)return(list("weight_bins"=drop_na(group_weights)))}```### Potential Journals:1. Global Change Biology2. ICES Journal of Marine Science# AbstractThe Northwest Atlantic Ocean one of the fastest warming locations on Earth. Research on the impacts of this rapid warming has primarily focused on high-profile and/or upper trophic level species. The temperature size rule (TSR) suggests that ectotherms living at elevated temperatures will reach smaller body-sizes at comparable development stages. However, it is unclear whether that relationship can be mitigated against at an individual level through adaptive behaviors in an open ocean environment. Here we've investigated ecosystem wide impacts using size spectrum analysis to track changes in the community size structure. In cases where community responses fail at adapting to elevated temperatures, we anticipated a steepening of the body-size spectrum slope relating to a reduction in larger sized individuals and an increased prevalence of smaller sized individuals. Using data from fisheries independent surveys we estimated size spectrum indices for four regions along the US NE continental shelf. Correlation/regression analyses were then performed to then assess the degree to which these changes were in alignment to hypothesized bottom-up and top-down disturbances. At the regional scale, we found that community size structure changes (spectra slope) were the largest in the Northern regions, in the Gulf of Maine and Georges Bank. These areas are home to the coldest temperatures and the largest proportions of groundfish species in the community. Spectrum slope declines were most pronounced in the 80's and 90's, before to the more-recent elevated temperatures. This suggests that external forces besides temperature drove the initial declines of larger-sized individuals within the communities, before elevated temperatures began to influence the ecosystem. Correlation analyses reveal that while fisheries landings are strongly correlated with these declines, bottom-up factors of zooplankton community metrics, Gulf Stream Index, and SST anomalies are also important. While the primary pressure of fisheries exploitation has declined dramatically over time, the recovery of larger-sized individuals has not been seen. That kind of recovery will likely depend on the elevated temperatures seen over the last decade.# Introduction## Temperature & EcologyTemperature plays a critical role on biological life impacting many of the chemical reactions that underpin basic physiological function. Temperature has direct and indirect impacts on critical biological functions including the acquisition of biomass through feeding, the rates of biomass loss through metabolism, and the rates at which individuals mature and develop. Because of these relationships, most species have evolved thermal preferences around which these chemical reactions are optimized. Species that are unable to maintain their thermal preferences internally must be able to follow their thermal preference in the environment through locomotion or adapt to less-favorable conditions through changes in behavior or risk metabolic costs in failing to do so. In an era of anthropogenic climate change, there is an expectation that many species will be displaced from historic habitats in their efforts to follow their thermal preferences. Recent research in marine environments has shown evidence of this as species are now shifting to higher latitudes and to deeper depths in the pursuit of more favorable conditions ( @kleisner2017 @pinsky2013 ). Others have suggested that temperature related impacts may not be seen through geographic distribution change, but through physiological changes, changes in seasonal phenology, or in dashed hopes of species recovery ( @daan2005 @miller2018 @pershing2018 @meyer-gutbrod2021 ).**Need to connect temperature to size here**The potential for elevated temperatures to impact the size structure of an ecosystem has implications for the ecosystem resilience in the face of climate change, as well as the blue economies & natural resource systems that rely upon their good health.## Size Spectrum TheorySize is a defining characteristic of all species and mediates many ecological interactions ( @brown2000 ). Body size impacts the mobility of an organism, its ability to evade predation, its ability to find success preying on other organisms, its efforts to locate and follow essential habitats, and the metabolic costs associated with each of these behaviors. Body size also mediates how vulnerable individuals are to features their immediate environment such as temperature and/or resource availability.Size structured environments are a fundamental organizational pattern emerging from these relationships that has been heavily researched. In the context of a strongly size-structured ecosystem, growth and maturity changes alter fitness and ultimately determine whether a species is successful in the given environment. Ecological theory is rich with models relating how energy transfers from smaller prey species to larger predatory trophic levels, the allocation of energy for growth, and the trade offs of allocating energy towards reproduction at the expense of growth. A persistent pattern in ecological systems globally is the relationship between body size and abundance.> \^ Add something about GSDR to just state what the relationship isSize spectrum and individual size distribution (ISD) models have gained increasing attention as an entry point to assessing ecosystem health and to detect systemwide disturbance ( @shin2005 ). These models avoid the need to explicitly articulate each predator-prey interactions and can be estimated from the commonly collected measures of abundance and size. A "size spectrum" describes the distribution of biomass or abundance as a function of individuals' mass or size on a log--log scale @guiet2016 . Size spectra are described by two terms, the size spectrum slope & intercept. These two terms convey the baseline productivity, and how energy flows through an ecosystem in the form of biomass. The spectrum intercept value captures the richness or the productivity at the base of the community and is strongly connected to the prevailing environmental conditions @boudreau1992. So much so, that spectrum shape is sometimes defined by its eutrophic or oligotrophic environmental conditions @rossberg2012.Size spectra condense the complexities of predator prey networks and their interactions into a handful of size-based indices (SBI). In doing so a community size spectrum captures the emergent properties of a system, while needing only size and abundance information that is commonly available. For this reason it has become increasingly relied upon as an indicator of ecosystem health in the push for ecosystem based fisheries management. Changes in slope have been associated with fishing exploitation, primarily through the targeted removal of larger individuals ( @bianchi2000; @shin2005 ). Numerical experiments have also tried to link changes in slope to environmental disturbances @guiet2016 . Biomass spectrum present predictable intercepts between ecosystems of similar productivity levels, but also of distinct temperatures @guiet2016.Direct quote from @guiet2016 , but nails the connection back to temp expectations:> Because it controls chemical reactions, temperature controls metabolic rates which underpin maintenance, growth or reproduction (Clarke and Johnston, 1999; Kooijman, 2010) as well as the functional responses to food density (Rall et al., 2012). Guiet et al. (2016)... In addition to the impact of temperature on communities' intercepts (heights), the impact of temperature on the speed of the energy flow within communities may affect other properties, such as their resilience to perturbations or the intensity of trophic cascades (Andersen and Pedersen, 2009).### Species Trends in the Northeast Atlantic ShelfThe continental shelf groundfish community in Northwest Atlantic has changed dramatically over the last century. Stocks that supported international fishing effort collapsed, and recovery efforts fell short of their objectives. Research on Georges Bank estimated that biomass more than halved in the 1960's (pre-dating federal monitoring efforts), and noted a species replacement of commercial groundfish target species by skate and dogfish @fogarty1998. Industrial fishing is inherently size-selective, with larger individuals selectively removed from the population. This has an immediate impact on the community size-distribution with additional impacts on the future population as well. Larger individuals have a greater impact on population recovery, capable of holding more (and often of higher quality) eggs. Size-based harvest in fisheries has been shown to create selective pressures that promote characteristics of early maturation at smaller sizes.### Temperature of the Gulf of Maine & NE ShelfIn addition to the ecological disturbance of industrial fisheries, the Northwest Atlantic is also one of the fastest warming locations in the global oceans. Sea surface temperatures in the Gulf of Maine since 1982 have been warming at rates faster than 96% of the world's oceans, with similar warming rates along the northwest Atlantic shelf @pershing2018. This persistent elevated temperature regime of the area is a result of several forces, a combination of shifting ocean currents and the unique bathymetry of the region. A Northward shift in the Gulf Stream directly increased the regional temperatures through increased transport of warm Gulf Stream water into areas like the Gulf of Maine. The Northward Gulf Stream shift is associated with a higher frequency of warm core rings, and the obstruction of cold-water Scotian Shelf current flow that would otherwise counter the influence of the Gulf Stream on the region's temperatures ( @gangopadhyay2019 @meyer-gutbrod2021 ). The combination of these oceanographic changes has led to a warmer continental shelf habitat.The rapid warming in the northwest Atlantic is a major factor in the redistribution of marine species along the US east coast. Species have responded by adjusting the timing and locations of their seasonal migrations and shifting their geographic ranges ( @nye2009 @staudinger2019 ). There is evidence that warming has hampered fisheries recoveries as well. Adding a metabolic tax to physiological pathways like growth and metabolism. Species like Black Sea Bass, Atlantic shortfin squid, and Blue crab have been high-profile examples of species expanding their ranges to follow their thermal preferences. While species like the American lobster have shown declines at their southern range near Long Island Sound, with much doubt whether they will recover under the present temperature trends. The recent regime shift in the physical oceanography has also shown to be a catalyst for biological shifts as well ( @meyer-gutbrod2021 ; @perretti2017 ).While these examples show that species can respond to changes in the physical environment around them through movement & behavior, research elsewhere suggests that physiological responses integrated across species will manifest as changes in community size structure.### PurposeWith the understanding that populations depend on the health of their ecosystems, there is a need to have community-wide metrics to effectively understand and manage marine resources. Size based indices are metrics that can be estimated from the information that has historically been available from long-term survey efforts. These indices have been shown to be sensitive to the impacts of fishing, but should also capture environmentally driven stress as well. We estimated size spectrum relationships as SBI's for the groundfish populations for each sub-region of the Northeast US continental shelf. In the case of the NW Atlantic sustained increases in temperature should have a physiological impact on the community size structure.This leads to our second hypothesis:> H2. Warming alters the community through the direct influence of temperature on metabolism, growth, and population productivity.# Methods## Groundfish Data```{r}#| label: load survdat_data# 1. Biological data used as inputwithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(catch_1g_labelled)) # rename and formatcatch_size_bins <- catch_1g_labelled# Get total number of speciesn_species <-length(unique(catch_size_bins$comname))```Fishery Independent data on was collected as part of the NEFSC's northeast trawl survey. This survey is conducted each year in the spring and in the fall, with sample locations determined following a stratified-random survey design with effort allocated in proportion to stratum area. Trawls are performed for a fixed duration at each station, reporting abundance at length and total biomass for all species caught. Trawl survey analyses incorporated data from both the Spring and Fall survey seasons, and for all years from 1970 to 2019. Correction factors were applied to total species abundance and aggregate species biomass to account for all changes in vessels, gear, and doors when appropriate as part of the survey program. However, abundance and biomass at length is not corrected for as part of the standard survey data protocol and needed to be estimated. To account for this, abundance at length for each species were adjusted to match the correction factors applied to total species abundance at each station, with allocations following the distributions of length caught at that station. Such that for each species: the sum of the resulting adjusted abundance numbers across each length is equal to the total abundance that was corrected for changes in vessels, gear, etc. To account for differences in sampling effort among survey strata, all corrected abundance-at-length data was then area-stratified.```{r}#| label: make region map# Get region polygons:# Get the paths to the shapefiles used to masktrawl_paths <- gmRi::get_timeseries_paths(region_group ="nmfs_trawl_regions", box_location ="cloudstorage")# Polygons for each regionall_poly <-read_sf(trawl_paths[["inuse_strata"]][["shape_path"]]) # Make the Mapggplot() +geom_sf(data = us_poly) +geom_sf(data = canada) +geom_sf(data = all_poly, aes(fill = survey_area), alpha =0.8) +coord_sf(xlim =c(-76.4, -64.4), ylim =c(35, 45.5), expand = F) +scale_fill_gmri() +theme_bw() +map_theme(legend.position ="bottom", legend.background =element_rect(fill ="white")) +guides(fill =guide_legend(title ="", nrow =2)) +labs(title ="Groundfish Survey Regions")```Data from the survey was grouped by strata to form geographically meaningful sub-regions: Gulf of Maine, Georges Bank, Southern New England, Mid-Atlantic Bight. For each region, we developed several time series indicators:1. Community Composition metrics (abundance and biomass by functional group, with body-size contributions)2. Mean size of the aggregate community and key functional groups3. Slope and intercept of the size spectrum### Community CompositionFunctional groups were assigned to each species based on life history and geography using the definitions of @hare2010. Functional groups included were coastal, diadromous, elasmobranch, groundfish, pelagic, and reef species. Stratified annual abundance and biomass totals were calculated for each functional group and each region with labels for increasing body-size (biomass kg) groups.```{r}#| label: biomass-group-summaries# The handmade weight and length bins# # Do some grouping by year survey area and functional group# fgroup_area <- get_group_summaries(data_in = catch_size_bins, Year, survey_area, hare_group)# fgroup_all <- get_group_summaries(data_in = catch_size_bins, Year, hare_group)# fgroup_all <- fgroup_all %>% map(~mutate(.x, survey_area = "Northeast Shelf, Full"))# Actual loog10 bins# Do some grouping by year survey area and functional groupfgroup_area <-l10_bin_summaries(data_in = catch_size_bins, Year, survey_area, hare_group)fgroup_all <-l10_bin_summaries(data_in = catch_size_bins, Year, hare_group)fgroup_all <- fgroup_all %>%map(~mutate(.x, survey_area ="Northeast Shelf, Full"))# Join the overall values back infgroup_area <-map2(fgroup_area, fgroup_all, bind_rows)```### Body Size TrendsAnnual abundance weighted body length and body weight within each region and for each functional group were also estimated using the numbers at length and the estimated biomass at length information. Data for body size trends were not truncated and include all sizes caught.```{r}#| label: body size data# Load the average body size datawithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(mean_sizes_ss_groups)) # Grouped on year and regionregional_size <- mean_sizes_ss_groups %>%filter(`group ID`=="single years * region") %>%mutate(Year =as.numeric(Year),survey_area =factor(survey_area, levels = area_levels))# Grouped on year, region, & functional groupfunctional_group_size <- mean_sizes_ss_groups %>%filter(`group ID`=="single years * region * functional group") %>%mutate(Year =as.numeric(Year),survey_area =factor(survey_area, levels = area_levels))```### Size Spectrum AnalysesCommunity size spectra were estimated using abundance-at-length data from `r n_species` species. These species were selected based on the availability of published weight-at-length relationships (Wigley et al. 2003) and represented 98.98% of the total biomass caught in the survey. Published length-weight relationships were used to convert abundance at length data into their corresponding biomass at length (kg). These values were then used to get totals for stratified weight-at-length, in complement to the corrected abundances-at-length data which had been area stratified. These area-stratified biomass at length totals were then used for fitting each regional biomass size spectra and exponents of individual size distributions. Data for these analyses was truncated by body masses of a minimum 1g and a maximum 10kg to account for poor gear selectivity at the smallest and largest size ranges and to ensure that all size bins were sampled in all years for the body mass spectrum fitting.#### Body Mass SpectraNormalized biomass spectra were estimated following the LBNbiom methodology as described in @edwards2017. When fitting the normalized biomass size spectra, stratified biomass at length data was binned into equal spaced intervals on a (1) on a $log_{10}$ scale, with bodymass totaled across all species. To normalize the spectra, the stratified abundances within each bin was then divided by the bin-width to account for the increasing bin-widths, a consequence of the log scale. Normalized bodymass spectra were fit for each year and for each region independently, and for each year across all strata, using ordinary least squares (OLS) regressions for stratified abundance (normalized) by body-size bins.```{r}#| label: load the size spectrum indiceswithr::with_dir(rprojroot::find_root('_targets.R'), tar_load(size_spectrum_indices)) # Grab SS Groups we care aboutregion_indices <- size_spectrum_indices %>%filter(`group ID`=="single years * region") %>%mutate(yr =as.numeric(as.character(Year)),survey_area =factor(survey_area, levels = area_levels),sig_fit =ifelse(l10_sig_strat <0.05, "Significant", "Non-Significant"))```#### Individual Size DistributionsIndividuals in the catch data are measured to the nearest cm, with smaller individuals measured to the nearest millimeter. Because individual biomass is estimated from those length measurements, there is for each increment a range of possible bodymass that spans between the cm & mm measurement increments. The relationships between length and bodymass in fishes is exponential and taxon specific, so biases resulting from using only the lower or upper end of those ranges is different for each taxon and increases for larger taxa. To account for this, we used the extended likelihood method (MLEbin) of @edwards2020 . This method estimates the exponent of size spectra (b) for a bounded power law relationship between individual abundances and their estimated biomass. With the MLEbin method the ranges of bodymass associated with each length measurement is accounted for to reduce bias. Exponents of the individual size distribution were estimated for each region and for every year. A minimum body mass of 1g was used for the lower bound with 10kg individual bodymass used as an upper bound for the ISD relationship. The biomass within these limits represents 97.83% of all estimated biomass used in these analyses.## Temperature DataGlobal Sea surface temperature data was obtained via NOAA's optimally interpolated SST analysis (OISSTv2), providing daily temperature values at a 0.25° latitude x 0.25° longitude resolution (Reynolds et al. 2007). A daily climatology for every 0.25° pixel in the global data set was created using average daily temperatures spanning the period of 1982-2011. Daily anomalies were then computed as the difference between observed temperatures and the daily climatological average. OISSTv2 data used in these analyses were provided by the NOAA PSL, Boulder, Colorado, USA from their website at https://psl.noaa.gov.Temperature data was regionally averaged to match the survey regions from the age-at-length data. SST anomalies were averaged by year for each region and over the entire sampling region to produce daily time series. These time series were then processed into annual timeseries of surface temperatures and anomalies. All region-averaging was done with area-weighting of the latitude/longitude grid cells to account for differences in cell-size in of the OISSTv2 data.```{r}#| label: load temperature data# Load the regional temperatures from {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(regional_oisst))# Get regional averagestemp_regimes <- regional_oisst %>%filter(yr >1981) %>%mutate(regime =ifelse(sst_anom >0, "hot", "cold"),survey_area =ifelse(survey_area =="all", "Northeast Shelf, Full", survey_area),survey_area =factor(survey_area, levels = area_levels)) ```## Spectra DriversThe impact of external factors on the changes in size spectra was correlated against several hypothesized driving forces related to both environmental regimes and anthropogenic disturbances. Potential environmental drivers include sea surface temperature anomalies, Gulf Stream Index (GSI), and two zooplankton community indices derived from the Gulf of Maine continuous plankton recorder (CPR) program. Hypothesized anthropogenic drivers include state and federal fisheries landings from the Greater Atlantic Regional Fisheries Office (GARFO), divided by reporting zones into aggregate regions to closely align with the survey areas we defined for the size spectra analyses.{fig-align="center" width="473"}```{r}#| label: spectra-driver-data#| eval: true#### 1. GARFO Landings ##### Load the GARFO landings data:res_path <- gmRi::cs_path("res")# Landings of finfish* sheet 5landings <-read_xlsx(path =str_c(res_path, "GARFO_landings/KMills_landings by area 1964-2021_JUN 2022.xlsx"), sheet =5) %>%rename_all(tolower)# Shapefiles for the fisheries stat zonesstat_zones <-read_sf(str_c(res_path, "Shapefiles/Statistical_Areas/Statistical_Areas_2010_withNames.shp"))# Make a list of zones to roughly match the survey areas:fish_zones <-list("Gulf of Maine"=c(511:515, 464, 465),"Georges Bank"=c(521, 522, 525, 561, 562),"Southern New England"=c(611, 612, 613, 616, 526, 537, 538, 539),"Mid-Atlantic Bight"=c(614:615, 621, 622, 625, 626, 631, 632))# Trim out what we don't need and label themstat_zones <- stat_zones %>%mutate(survey_area =case_when( Id %in% fish_zones$"Gulf of Maine"~"Gulf of Maine", Id %in% fish_zones$"Georges Bank"~"Georges Bank", Id %in% fish_zones$"Southern New England"~"Southern New England", Id %in% fish_zones$"Mid-Atlantic Bight"~"Mid-Atlantic Bight")) %>%filter(survey_area %in%c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight"))# # Map it out? - takes forever from Chesapeake bay# garfo_map <- ggplot(stat_zones) +# geom_sf(aes(fill = survey_area), alpha = 0.8) +# geom_sf(data = us_poly) +# geom_sf(data = canada) +# coord_sf(xlim = c(-76.4, -64.4), ylim = c(35, 45.5), expand = F) +# scale_fill_gmri() +# theme_bw() +# map_theme(legend.position = "bottom",# legend.background = element_rect(fill = "white")) +# guides(fill = guide_legend(title = "", nrow = 2)) +# labs(title = "GARFO Region Assignments")# ggsave(garfo_map, filename = here::here("R/nmfs_size_spectra/images/garfo_regions.png"), bg = "white")# Add the labels into the landings data and remove what we don't need there:landings <- landings %>%mutate(survey_area =case_when(`stat area`%in% fish_zones$"Gulf of Maine"~"Gulf of Maine",`stat area`%in% fish_zones$"Georges Bank"~"Georges Bank",`stat area`%in% fish_zones$"Southern New England"~"Southern New England",`stat area`%in% fish_zones$"Mid-Atlantic Bight"~"Mid-Atlantic Bight")) %>%filter(survey_area %in%c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight")) %>%mutate(survey_area =factor(survey_area, area_levels_long))# Get Summarieslandings_summ <- landings %>%drop_na() %>%rename("weight_lb"=`landed lbs`,"live_lb"=`live lbs`) %>%drop_na() %>%group_by(year, survey_area) %>%summarise( across(.cols =c(value, weight_lb, live_lb), .fns =list(mean = mean, total = sum), .names ="{.fn}_{.col}"), .groups ="drop") # Scale the landings to create an index by arealandings_summ <- landings_summ %>%group_by(survey_area) %>%mutate(total_wt_z = base::scale(total_weight_lb)) %>%ungroup()#### 2. Climate Drivers ##### From ecodata:#### NAO & GSI & WCR# North Atlantic Oscillationnao <- ecodata::nao %>%mutate(year = Time,Time =as.Date(str_c(Time, "-01-01")))# Gulf Stream Indexgsi <- ecodata::gsi %>%mutate(Time =str_pad(as.character(Time), width =7, side ="right", pad ="0"),year =as.numeric(str_sub(Time, 1, 4)),month =str_sub(Time, -2, -1),Time =as.Date(str_c(year, month, "01", sep ="-")))# # Warm Core Rings# wcr<- ecodata::wcr %>%# mutate(# year = Time,# Time = as.Date(str_c(Time, "-01-01")))#### 3. CPR Indices ##### From CPR Analysis:# Uses 7 focal taxa, 2 calanus stagescpr_indices <-read_csv("~/Documents/Repositories/continuous_plankton_recorder/results_data/cpr_focal_pca_timeseries_period_1961-2017.csv")# Do some Reshaping for the Plotcpr_indices <- cpr_indices %>%select(c(year, 2:3)) %>%rename(`Principal Component 1`="First Mode",`Principal Component 2`="Second Mode") %>%pivot_longer(names_to ="pca_mode", values_to ="pca_loading", cols =c(2:3))# # Zooplankton densities from ECOMON?# # tidy the calanus stage abundance data# cal <- ecodata::calanus_stage %>% # rename(year = Time,# var = Var,# area = EPU) %>% # group_by(var, area) %>% # mutate(value = scale(Value)) %>% # ungroup() %>% # select(-Value)``````{r}#| label: all-drivers-table#### Build a Common form for all Drivers: ##### Format: year, area, var, value# table to join for swapping shorthand for long-hand namesarea_df <-data.frame(area =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "All"),survey_area =c("SS", "GoM", "GB", "SNE", "MAB", "Full Region"))#--------------------# 1.Climate Features# Join the climate modes togetherclim_drivers <-bind_rows(list(#nao, # Only Goes to 2018, also not strongly correlated gsi#, #wcr # Too related to GSI )) %>%filter(EPU =="All")# Put climate drivers on annual scheduleclim_idx <- clim_drivers %>%group_by(year, area = EPU, var = Var) %>%summarise(value =mean(Value, na.rm = T),.groups ="drop")#--------------------# 2. GARFO landingslandings_idx <- landings_summ %>%group_by(year, area = survey_area) %>%summarise(value =mean(total_wt_z, na.rm = T),.groups ="drop") %>%mutate(var ="landings")#--------------------# 3. CPR Zooplanktoncpr_idx <- cpr_indices %>%mutate(area ="Gulf of Maine",var =str_c("cpr_pca_", str_sub(pca_mode, -1, -1))) %>%group_by(year, area, var) %>%summarise(value =mean(pca_loading),.groups ="drop")#--------------------# 4. Temperaturesst_idx <- temp_regimes %>%left_join(area_df, by ="survey_area") %>%select(year = yr, area, value = sst_anom) %>%mutate(var ="sst_anom")#--------------------# 5. Size Spectrum Slopes & Interceptsss_idx <- region_indices %>%select(year = Year, survey_area, spectra_slope = l10_slope_strat, spectra_int = l10_int_strat,isd_slope = b) %>%pivot_longer(cols =c(spectra_slope, spectra_int, isd_slope), names_to ="var", values_to ="value") %>%mutate(year =as.numeric(year)) %>%left_join(area_df, by ="survey_area") %>%select(-survey_area)#--------------------# # ECODATA calanus# Not strongly correlated, also a big mess to handle# Less likely to be important relative to the broader zooplankton landscape# # tidy the calanus stage abundance data# # Scale the abundances within each area# cal_idx <- ecodata::calanus_stage %>% # mutate(EPU = ifelse(EPU == "GOM", "GoM", EPU)) %>% # group_by(Var, EPU) %>% # mutate(value = as.numeric(scale(Value))) %>% # ungroup() %>% # select(-Value, -Units) %>% # setNames(c("year", "var", "survey_area", "value")) %>% # left_join(area_df, by = "survey_area") %>% # select(-survey_area)#--------------------# All Metrics Togetherall_drivers <-bind_rows(list( clim_idx, landings_idx, cpr_idx, sst_idx, ss_idx ))``````{r}#| label: set-driver-periods#### Working Here ##### Reshape Wide# Scale the columnsdrivers_scaled <- all_drivers %>%mutate(id =str_replace_all(str_c(area, "_", var), " ", "_")) %>%select(-c(area, var)) %>%pivot_wider(names_from = id, values_from = value) %>%column_to_rownames(var ="year") %>%scale() %>%as.data.frame()# Pull different years into a listdriver_periods <-list("regime_1"=c(1954:2005),"regime_2"=c(2006:2020),"full"=c(1954:2020))# Split them out into a list that contains data for each perioddriver_period_l <-map(driver_periods, function(yrs){ period_matrix <- drivers_scaled[rownames(drivers_scaled) %in% yrs, ] %>%drop_na()})``````{r}#| label: debug-corrplot-organization#| eval: false# Recombine SST and Landings to be match the different areas # Get significance values# flag them on correlation plot# Grab scatterplots of the relationships too# Step 1.# Do some reshaping to create a matrix for correlations:d_test <- driver_period_l$regime_1# Step 1.# Get the correlation matrix (data is pre-scaled)corr_matrix <-rcorr(as.matrix(d_test))rho <- corr_matrix$r# Step 2. Significancesig_mat <- corr_matrix$P# Step 3. # Make correlation matrix & significance a dataframe for GGplot# Correlation Matrixgg_rho <-as.data.frame(rho) %>%rownames_to_column(var ="xvar") %>%pivot_longer(cols =-xvar, names_to ="yvar", values_to ="r2")# Significance Dataframegg_sig <-as.data.frame(sig_mat) %>%rownames_to_column(var ="xvar") %>%pivot_longer(cols =-xvar, names_to ="yvar", values_to ="pval")# Step 4. # Tidy it up, enhance the labels for plotting# label the drivers and the responses specifically so we can group them# in an orderly waygg_rho_tidy <-left_join(gg_rho, gg_sig) %>%mutate(# A. Slope and intercept are the only size spectra features# and alsoresponse_var =ifelse(str_detect(xvar, "slope|int"), T, F),# B.These are the drivers, use string detect to pull them outdriver_var =ifelse(str_detect(yvar, "adt|c5|osc|pca|landings|index|anom"), T, F),# C. Label Top down and Bottom Up Driversdriver_type =case_when(str_detect(yvar, "landings") ~"Top-Down",str_detect(yvar, "adt|c5") ~"Calanus",str_detect(yvar, "pca|stream|anom|osci") ~"Bottom-Up"),# D. These are the areas associated with the Spectra Featuresxvar_area =case_when(str_detect(xvar, "All") ~"All",str_detect(xvar, "Georges") ~"Georges Bank",str_detect(xvar, "Gulf") ~"Gulf of Maine",str_detect(xvar, "Southern") ~"Southern New England",str_detect(xvar, "Mid-Atlantic") ~"Mid-Atlantic Bight"),# E. Set Factor Levels for plotxvar_area =factor( xvar_area, levels =c("All", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")),# F. These are the areas of the driversyvar_area =case_when(str_detect(yvar, "All") ~"All",str_detect(yvar, "Georges") ~"Georges Bank",str_detect(yvar, "Gulf") ~"Gulf of Maine",str_detect(yvar, "Southern") ~"Southern New England",str_detect(yvar, "Mid-Atlantic") ~"Mid-Atlantic Bight"),yvar_area =factor(yvar_area, levels =c("All", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")),# Shorthand for the Spectra Parametersx_short =case_when(str_detect(xvar, "isd") ~"ISD Exponent",str_detect(xvar, "slope") ~"Spectra Slope",str_detect(xvar, "int") ~"Spectra Intercept"))# Match SST and Landings to their areas:# Rename landings and SST so they occupy the same plot space# Do some filtering and string formatting:only_driver_response <- gg_rho_tidy %>%filter( response_var, driver_var) %>%mutate(area_match = xvar_area == yvar_area) %>%filter(area_match |!str_detect(yvar, "sst|landings")) %>%mutate(yvar =ifelse(str_detect(yvar, "landings"), "Commercial_Landings", yvar),yvar =ifelse(str_detect(yvar, "sst"), "SST_Anomalies", yvar),yvar =ifelse(str_detect(yvar, "gulf_stream_index"), "Gulf Stream Index", yvar),yvar =str_replace_all(yvar, "_", " "),yvar =factor(yvar, levels =c("Gulf of Maine cpr pca 1","Gulf of Maine cpr pca 2","Gulf Stream Index","All north atlantic oscillation","SST Anomalies","Commercial Landings")) ) #------- # Do some flagging for significanceonly_driver_response <- only_driver_response %>%mutate(signif =ifelse(pval <0.05, "*", ""))# Check Correlation Plotonly_driver_response %>%ggplot(aes(yvar, x_short, fill = r2)) +geom_tile() +geom_text(aes(label = signif)) +facet_grid(xvar_area~driver_type, scales ="free_x", space ="free") +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +coord_cartesian(clip ="off") +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position =c(1.205, -.25),legend.title =element_text(vjust =0.75)) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5)) +labs(x ="Hypothesized Driver", y ="Bodymass Spectra Coefficient",fill ="Correlation Coefficient",title ="Correlation of Size Spectra to Hypothesized Drivers | All Years",subtitle ="Correlation coefficients shown display same-year correlations.",caption =md("* Symbol signals a significant correlation at an alpha = 0.05"))# Fine... Correlation Matrix Plots.d_test %>%select(matches("george|cpr|index", ignore.case =TRUE)) %>%rename_all(.funs = str_replace_all, pattern ="_", replacement =" ") %>%ggplot(aes(x = .panel_y, y = .panel_x)) +geom_point(alpha =0.4, shape =16) +geom_smooth(aes(colour =NULL, fill =NULL), method ="lm") +facet_matrix(cols =vars(-matches("int|slope")),rows =vars(matches("slope|int"))) +labs(x ="Spectra Coefficient (z)",y ="Driver (z)")``````{r}#| label: period-correlations-new# Run Correlations for Each period and for the entire periodperiod_correlations <-map(driver_period_l, function(drivers){# Step 1.# Get the correlation matrix (data is pre-scaled) corr_matrix <-rcorr(as.matrix(drivers)) rho <- corr_matrix$r# Step 2. Significance sig_mat <- corr_matrix$P# Step 3. # Make correlation matrix & significance a dataframe for GGplot gg_rho <-as.data.frame(rho) %>%rownames_to_column(var ="xvar") %>%pivot_longer(cols =-xvar, names_to ="yvar", values_to ="r2") gg_sig <-as.data.frame(sig_mat) %>%rownames_to_column(var ="xvar") %>%pivot_longer(cols =-xvar, names_to ="yvar", values_to ="pval")# Step 4. # Tidy it up, enhance the labels for plotting# label the drivers and the responses specifically so we can group them# in an orderly way gg_rho_tidy <-left_join(gg_rho, gg_sig, by =c("xvar", "yvar")) %>%mutate(# A. Slope and intercept are the only size spectra featuresresponse_var =ifelse(str_detect(xvar, "slope|int"), T, F),# B.These are the drivers, use string detect to pull them outdriver_var =ifelse(str_detect(yvar, "adt|c5|osc|pca|landings|index|anom"), T, F),driver_type =case_when(str_detect(yvar, "landings") ~"Top-Down",str_detect(yvar, "adt|c5") ~"Calanus",str_detect(yvar, "pca|stream|anom|osci") ~"Bottom-Up"),# C. These are the areas associated with the Spectra Featuresxvar_area =case_when(str_detect(xvar, "All") ~"All",str_detect(xvar, "Georges") ~"Georges Bank",str_detect(xvar, "Gulf") ~"Gulf of Maine",str_detect(xvar, "Southern") ~"Southern New England",str_detect(xvar, "Mid-Atlantic") ~"Mid-Atlantic Bight"),xvar_area =factor( xvar_area, levels =c("All", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")),# D. These are the areas of the driversyvar_area =case_when(str_detect(yvar, "All") ~"All",str_detect(yvar, "Georges") ~"Georges Bank",str_detect(yvar, "Gulf") ~"Gulf of Maine",str_detect(yvar, "Southern") ~"Southern New England",str_detect(yvar, "Mid-Atlantic") ~"Mid-Atlantic Bight"),yvar_area =factor( yvar_area, levels =c("All", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")),# Shorthand for the Spectra Parametersx_short =case_when(str_detect(xvar, "isd") ~"ISD Exponent",str_detect(xvar, "slope") ~"Spectra Slope",str_detect(xvar, "int") ~"Spectra Intercept"),x_short =factor(x_short, levels =c("Spectra Intercept","Spectra Slope","ISD Exponent")))# Match SST and Landings to their areas:# Rename landings and SST so they occupy the same plot space# Do some filtering and string formatting: only_driver_response <- gg_rho_tidy %>%filter( response_var, driver_var) %>%mutate(area_match = xvar_area == yvar_area) %>%filter(area_match |!str_detect(yvar, "sst|landings")) %>%mutate(yvar =ifelse(str_detect(yvar, "landings"), "Commercial_Landings", yvar),yvar =ifelse(str_detect(yvar, "sst"), "SST_Anomalies", yvar),yvar =ifelse(str_detect(yvar, "gulf_stream_index"), "Gulf Stream Index", yvar),yvar =str_replace_all(yvar, "_", " "),yvar =factor(yvar, levels =c("Gulf of Maine cpr pca 1","Gulf of Maine cpr pca 2","Gulf Stream Index","All north atlantic oscillation","SST Anomalies","Commercial Landings")),# Flag whether the relationship wa significant or notsignif =ifelse(pval <0.05, "*", "") ) # Spit out all the pieces in an organized list correlation_details <-list("drivers_matrix"= drivers,"rho"= rho,"pvals"= sig_mat,"ggrho"= gg_rho,"ggrho_tidy"= gg_rho_tidy,"only_driver_response"= only_driver_response )})```# Results### Abundance DistributionAbundance across all body sizes remained relatively stable from the 1970's before rising in the northern regions around 1990 beginning in the Gulf of Maine. Around this time abundances increased through the mid 2010's. Further south in Georges Bank, abundances remained flat until the 2010's, when overall abundance roughly tripled, only to fall back to previous amounts by the end of the century. Southern New England saw a less dramatic rise and fall that began just before 2010, again falling back to earlier levels by the end of the century. The Mid Atlantic Bight had relatively consistent abundances throughout, with no major periods of abundance growth or decline, but with larger inter-annual variability.```{r}#| label: abundance distributions#| fig.width: 8#| fig-height: 8# panel plot the biomass by body sizefgroup_area$weight_bins <- fgroup_area$weight_bins %>%mutate(strat_abund_mill = wtbin_strat_abund /1e6,strat_lwbio_mill = wtbin_strat_lw_bio /1e6,survey_area =factor(survey_area,levels = area_levels))# Just abundance, not by groupsstrat_abundance_plot <- fgroup_area$weight_bins %>%group_by(Year, survey_area, hare_group) %>%summarise(abund_millions =sum(strat_abund_mill, na.rm = T),.groups ="drop") %>%ggplot(aes(Year, abund_millions, fill = hare_group)) +geom_area(color ="gray30", size =0.1) +scale_fill_gmri() +scale_x_continuous(expand =expansion(add =c(0,0))) +facet_wrap(~survey_area, ncol =1, scales ="free") +labs(y ="Abundance (millions)", title ="Estimated Abundance from\nArea-Stratified Catch Rates",x ="",fill ="Functional Group") +theme(legend.position ="bottom",panel.background =element_rect(colour ="black", size=2, fill =NA))# Just Biomass, by area and functional groupstrat_biomass_plot <- fgroup_area$weight_bins %>%group_by(Year, survey_area, hare_group) %>%summarise(abund_millions =sum(strat_lwbio_mill, na.rm = T),.groups ="drop") %>%ggplot(aes(Year, abund_millions, fill = hare_group)) +geom_area(color ="gray30", size =0.1) +scale_fill_gmri() +scale_x_continuous(expand =expansion(add =c(0,0))) +facet_wrap(~survey_area, ncol =1, scales ="free") +labs(y ="Biomass (million kg)", title ="Estimated Biomass from\nArea-Stratified Catch Rates",x ="",fill ="Functional Group") +theme(legend.position ="bottom",panel.background =element_rect(colour ="black", size=2, fill =NA))# Biomass and Abundance Together(strat_abundance_plot | strat_biomass_plot) +plot_layout(guides ="collect") &theme(legend.position ="bottom")```Abundance gains observed in Georges Bank and Gulf of Maine were primarily from groundfish species, with additional growth from demersal species in the Gulf of Maine. Increases in abundance across all areas was mostly confined to individuals weighing less than .5kg. With some years driven in large-part by exceptional year-classes in a handful of species like haddock in Georges Bank. The observed abundance volatility in Southern New England and the Mid-Atlantic Bight conversely was largely the result of changes in abundance in pelagic species, whose abundance varied by several times the magnitude that of the other functional groups.```{r}#| label: abundance-by-group#| fig.height: 8#| fig.width: 10#| eval: true# Plot abundance by region and Groupfgroup_area$weight_bins %>%# ggplot(aes(Year, strat_abund_mill, fill = weight_bin)) +ggplot(aes(Year, strat_abund_mill, fill =factor(left_lim))) +geom_area(color ="transparent") +facet_grid(hare_group~survey_area, scales ="free") +scale_fill_gmri() +scale_y_continuous(labels =comma_format()) +labs(y ="Stratified Abundance (millions)", title ="Abundance Allocation by Weight",x ="", fill ="Minimum Individual Weight 10^x (g)") +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1))```### Biomass DistributionOverall biomass was highest in the two northern regions, the Gulf of Maine and Georges Bank. Roughly half of the biomass sampled in these regions can be attributed to groundfish/demersal species, with the second largest contributions coming from elasmobranchs. Groundfish biomass, larger individuals \>2kg in particular, declined during the 70's and 80's in these regions, never truly recovering. Beginning in the 2000's there were signs that groundfish abundances were increasing as evidenced by increasing numbers of smaller individuals, however in both regions this trend appears to have reversed by the mid 2010's. Elasmobranch biomass increased steadily throughout the survey time period across all regions, with the exception of southern New England. This area showed large 5-10 year swings in biomass, but no clear long-term trend. Larger elasmobranch were rare in all regions except for a period spanning the late 70's through the early 90's isolated to Georges Bank. Demersal species biomass was highest in the Gulf of Maine, dwarfing their contributions in other regions. Their biomass declined in the 70's, was flat until the late 90's, remaining relatively high until declining in the late 2010's. Pelagic species biomass was low in all regions, and is unlikely to be representative of true biomass trends due to gear selectivity.```{r}#| label: bodymass distributions#| fig.width: 10#| fig-height: 8# panel plot the biomass by body sizefgroup_area$weight_bins %>%complete(Year, hare_group, survey_area, left_lim) %>%# ggplot(aes(Year, strat_abund_mill, fill = weight_bin)) +ggplot(aes(Year, strat_abund_mill, fill =factor(left_lim))) +geom_area(color ="transparent") +facet_grid(hare_group~survey_area, scales ="free") +scale_fill_gmri() +scale_y_continuous(labels =comma_format()) +labs(y ="Stratified Total Biomass (million kg)", title ="Biomass Allocation by Weight",x ="", fill ="Minimum Individual Weight 10^x (g)") +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1))```### Regional Variation in Species CompositionThere was a distinct difference between Northern and Southern regions in the way biomass was distributed among the different functional groups. The primary contributors to overall biomass in the southern regions (southern New England & mid-Atlantic bight) was the elasmobranch community. While the northern regions (Gulf of Maine & Georges Bank) each had similar quantities of elasmobranch biomasses, there was also a comparable contribution of groundfish and in the Gulf of Maine there was a major component of demersal species as well.### Body Size TrendsThe average fish size in the Gulf of Maine (length and weight) declined the greatest of all regions over our study period. The average individual length was greatest in the 1970's in the 35-40cm range, falling to 28-33cm over the last decade. Body-weight fell dramatically in the 1980's, from around .75kg in the 1970's to .25-.30kg, roughly a third of what it had been. Georges Bank body sizes also declined during the study period, but less dramatically. Both of these Northern regions had brief period in the late 2000's where average length and weight rose, before falling again in the 2010's. The MAB region was the only region to see a long-term increase in both length and weight during the study period. SNE saw no long-term change in length, and a minor decline in average body-weight.::: panel-tabset#### All Species```{r}#| label: regional-size-trends#| fig.height: 8#| fig.weight: 8# Re-factorregional_size <- regional_size %>%mutate(survey_area =factor(survey_area, area_levels),area_copy = survey_area)# Length plotavg_len_p <- regional_size %>%ggplot(aes(Year, mean_len_cm)) +geom_line(data =select(regional_size, -survey_area),aes(group = area_copy),alpha =0.2, size =0.5) +geom_line(aes(group = survey_area#,#color = area_copy ), size =1,show.legend = F) +scale_color_gmri() +facet_wrap(~survey_area, ncol =1) +labs(title ="Average Length",y ="Length (cm)",color ="Region")# Weight plotavg_wt_p <- regional_size %>%ggplot(aes(Year, mean_wt_kg)) +geom_line(data =select(regional_size, -survey_area),aes(group = area_copy),alpha =0.2, size =0.5) +geom_line(aes(group = survey_area#,#color = area_copy ), size =1,show.legend = F) +scale_color_gmri() +facet_wrap(~survey_area, ncol =1) +labs(title ="Average Weight",y ="Weight (kg)",color ="Region")# Plot side by side(avg_len_p | avg_wt_p) ```#### Functional Group Length```{r}#| label: functional-group-length-trends#| fig.width: 10#| fig-height: 8# Refactor areasfunctional_group_size <- functional_group_size %>%mutate(survey_area =factor(survey_area, area_levels),area_copy = survey_area)functional_group_size %>%ggplot(aes(Year, mean_len_cm)) +geom_line(aes(group = survey_area#,#color = area_copy ), size =1,show.legend = F) +facet_grid(hare_group~survey_area, scales ="free_y") +labs(title ="Average Length",y ="Length (cm)",color ="Region")```#### Functional Group Bodymass```{r}#| label: functional-group-weight-trends#| fig.width: 10#| fig-height: 8# Average individual weightfunctional_group_size %>%ggplot(aes(Year, mean_wt_kg)) +geom_line(aes(group = survey_area), size =1,show.legend = F) +facet_grid(hare_group~survey_area, scales ="free_y") +labs(title ="Average Weight",y ="Weight (kg)",color ="Region")```:::### Regional Size SpectraBeginning in the 1970's, there was a clear difference in the relative positions of spectra parameters among the different regions. Gulf of Maine and Georges Bank showed the least steep spectra slopes in the earlier time periods with slopes around -1 & -1.1 respectively. The relatively flat slopes in these regions both steepened over time, settling near -1.3 (GoM) and -1.5 (GB). Gulf of Maine experienced much of its decline during the 1980's and 1990's. There was a brief reversal in this trend during the 2000's, but slopes continued to steepen by 2010 and remained steep through 2019. Georges Bank did not experience as rapid of a decline, but experienced a similar long-term steepening. In contrast to the northern regions, SNE and MAB had steeper slopes in the -1.2 to -1.5 territory. The long term pattern for SNE was one of increasing volatility, but not so much a decline. The spectra slope for the MAB was less volatile, but similarly maintained a relatively stable wander around -1.4. By the end of the study period all regions had slopes that were at or near a similar level.```{r}#| label: size spectra results#| fig.height: 8#| fig.weight: 8# Make a copy so we can gray out the linesregion_indices <-mutate(region_indices, area_copy = survey_area)# Plot the Regional Slopesslope_timeline <- region_indices %>%ggplot(aes(yr, l10_slope_strat, group = survey_area)) +geom_line(data =select(region_indices, -survey_area),aes(group = area_copy),alpha =0.2, size =0.5) +geom_line(aes(group = survey_area), size =1) +scale_color_gmri() +facet_wrap(~survey_area, ncol =1) +# theme_minimal() +labs(x ="Year", y ="Slope", title ="Normalized Biomass Spectra") # Plot the Regional Interceptsint_timeline <-ggplot(region_indices, aes(yr, l10_int_strat, group = survey_area)) +geom_line(data =select(region_indices, -survey_area),aes(group = area_copy),alpha =0.2, size =0.5) +geom_line(aes(group = survey_area), size =1) +facet_wrap(~survey_area, ncol =1) +# theme_minimal() +labs(x ="Year", y ="Intercept") # Assembleslope_timeline | int_timeline ```### Individual Size Distribution```{r}# Plot the Regional Slopesisd_timeline <- region_indices %>%ggplot(aes(yr, b, group = survey_area)) +geom_line(data =select(region_indices, -survey_area),aes(group = area_copy),alpha =0.2, size =0.5) +geom_line(aes(group = survey_area), size =1) +#geom_pointrange(aes(group = survey_area, ymin = confMin, ymax = confMax), size = 0.1) +scale_color_gmri() +facet_wrap(~survey_area, ncol =2) +labs(x ="Year", y ="Slope", title ="Normalized Biomass Spectra") # Assembleisd_timeline ```## Size Spectra Drivers### Driver Correlations**NOTE:** Correlation matrix is computed starting at the year where there are no NA values in any drivers. Currently with SST included that begins the matrix at 1982.::: panel-tabset#### Regime 1 (1982* - 2005)```{r}#| label: ggcorr-drivers-regime_1#| fig-height: 6#| fig-width: 8# Plot the convoluted thing:period_correlations$regime_1$only_driver_response %>%ggplot(aes(yvar, x_short, fill = r2)) +geom_tile() +geom_text(aes(label = signif)) +facet_grid(xvar_area~driver_type, scales ="free_x", space ="free") +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +coord_cartesian(clip ="off") +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position =c(1.205, -.25),legend.title =element_text(vjust =0.75),plot.caption =element_text(hjust =0)) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5)) +labs(x ="Hypothesized Driver", y ="Bodymass Spectra Coefficient",fill ="Correlation Coefficient",title ="Correlation of Size Spectra to Hypothesized Drivers | All Years",subtitle ="Correlation coefficients shown display same-year correlations.",caption =md("* Symbol signals a significant correlation at an alpha = 0.05"))```#### Regime 2 (2006-2019)```{r}#| label: ggcorr-drivers-regime_2#| fig-height: 6#| fig-width: 8# Plot the convoluted thing:period_correlations$regime_2$only_driver_response %>%ggplot(aes(yvar, x_short, fill = r2)) +geom_tile() +geom_text(aes(label = signif)) +facet_grid(xvar_area~driver_type, scales ="free_x", space ="free") +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +coord_cartesian(clip ="off") +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position =c(1.205, -.25),legend.title =element_text(vjust =0.75),plot.caption =element_text(hjust =0)) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5)) +labs(x ="Hypothesized Driver", y ="Bodymass Spectra Coefficient",fill ="Correlation Coefficient",title ="Correlation of Size Spectra to Hypothesized Drivers | Regime 2",subtitle ="Correlation coefficients shown display same-year correlations.",caption =md("* Symbol signals a significant correlation at an alpha = 0.05"))```#### Full Time Series```{r}#| label: ggcorr-drivers-full#| fig-height: 6#| fig-width: 8# Plot the convoluted thing:period_correlations$full$only_driver_response %>%ggplot(aes(yvar, x_short, fill = r2)) +geom_tile() +geom_text(aes(label = signif)) +facet_grid(xvar_area~driver_type, scales ="free_x", space ="free") +scale_fill_distiller(palette ="RdBu", limits =c(-1, 1), breaks =c(-1,0,1)) +coord_cartesian(clip ="off") +theme(axis.text.x =element_text(angle =45, hjust =1, vjust =1),strip.text.y =element_text(angle =0),legend.position =c(1.205, -.25),legend.title =element_text(vjust =0.75),plot.caption =element_text(hjust =0)) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5)) +labs(x ="Hypothesized Driver", y ="Bodymass Spectra Coefficient",fill ="Correlation Coefficient",title ="Correlation of Size Spectra to Hypothesized Drivers | All Years",subtitle ="Correlation coefficients shown display same-year correlations.",caption =md("* Symbol signals a significant correlation at an alpha = 0.05"))```:::# Discussion- Top-down and bottom up influences on both carrying capacity (intercept) and transfer efficiency (slope)Some of the major drivers suggested here operate on both, but to varying degrees. Here are some potential mechanisms:**Literature suggests:** - Intercept (a proxy for productivity and carrying capacity) is primarily determined by bottom up features like: nutrient availability, temperature- Slope (a measure of energy transfer efficiency and static biomass distribution) has been shown to be sensitive to the physical removal of species through fishing.**Temperature Mechanisms:** - Temperature's impact on growth via genetic plasticity impacts both the available biomass at the primary producer level (impacting ecosystem carrying capacity), as well as the Linf of larger species (recruitment rate). - Temperature also impacts the efficiency of energy (as biomass) being transferred between individuals via predator & prey interactions. More energy per-capita is expended in the form of increased metabolic rates and/or behavioral changes. This metabolic tax should steepen the spectrum slope by removing available energy at a system wide level. - Temperature impacts behavior via physiological impacts on metabolism and foraging rates as well as through the avoidance of temperature stresses.## Separating Complimentary Forces Impacting GrowthThe data used for this analysis was collected as part of a survey program which began out of concern that fisheries were already being over-harvested. Early estimates from scientists at that time suggested that by the 1970's total biomass of Georges Bank had already been halved, and elasmobranch species had begun to replace the over-harvested gadoid species [@fogarty1998]. Having such a large disturbance which pre-dates our time series is suggestive that the measured steepening of size spectrum slope we observed in this area and the adjacent Gulf of Maine are potentially the tail-ends of a longer and more severe ecosystem decline. While metrics of overall fishing pressure do not align exactly with trawl survey coverage, historical records and anecdotal evidence fro that time suggest that groundfish fishing pressure in these areas are a fraction of their what their impacts were in the 1960's and 1970's.## Forces Preventing the Recovery of Large IndividualsThis begs the question of why larger adult numbers never began to recover in these regions. Looking at abundance and biomass information from the survey there was evidence of strong recruitment among smaller individuals \< 1kg, but there has since not been a recovery in fishes larger than 1kg outside of the elasmobranchs. Work by [@pershing2015] suggested this prolonged recovery period may be due to a lack of accounting for temperature change in fisheries management. At the time of that research, the regional temperatures had only begun rising, and could have been considered at that time an acute stressor. Since that time the region has experienced nearly a decade of sustained above-average temperatures. There are signs that the success seen in recruitment and survival of even the smaller size classes is declining. While temperature change has been associated with changes in growth rates and size-at-age, so too have size-selective fishing practices, making it difficult to disentangle the importance of exploitation & temperature on the overall community size structure when body size integrates these two forces [@shackell2007].### Potential Drivers Timeseries:::: panel-tabset#### Surface Temperature Anomalies```{r}#| label: temperature-index#| fig.height: 8# Plot the timelinestemp_regimes %>%mutate(anom_direction =ifelse(sst_anom >0, "Positive", "Negative")) %>%ggplot( aes(yr, sst_anom)) +geom_col(aes(fill = anom_direction), size =0.75, alpha =0.8) +geom_hline(yintercept =0, size =0.5, lty =1) +scale_x_continuous(expand =expansion(add =c(0.25,0.25))) +scale_y_continuous(labels =number_format(suffix =" \u00b0C")) +scale_fill_gmri() +facet_wrap(~survey_area, ncol =1) +theme(legend.title =element_blank(),legend.background =element_rect(fill ="transparent"),legend.key =element_rect(fill ="transparent", color ="transparent")) +labs(title ="NW-Atlantic Regional SST Anomalies",x ="", y ="Surface Temperature Anomaly",caption ="Anomalies calculated using 1982-2011 reference period.") +guides(fill ="none",label ="none",color =guide_legend(keyheight =unit(0.5, "cm")))```#### GARFO Landings```{r}#| label: landings-index#| fig.height: 8#### Landings ####landings_summ %>%mutate(anom_direction =ifelse(total_wt_z >0, "Positive", "Negative")) %>%ggplot(aes(year, total_wt_z, group = survey_area)) +geom_col(aes(fill = anom_direction), size =0.75, alpha =0.8) +geom_hline(yintercept =0, size =0.5, lty =1) +scale_fill_gmri() +facet_wrap(~survey_area, ncol =1) +scale_x_continuous(breaks =seq(1950, 2030, by =10)) +labs(x ="Year", y ="Finfish Landings Index (z)",title ="GARFO Finfish Landings") +guides(fill ="none")```#### Climate Modes```{r}#| label: climate-modes#### Climate Drivers ####clim_drivers %>%group_by(year = lubridate::year(Time), Var, EPU) %>%summarise(Value =mean(Value, na.rm = T)) %>%mutate(anom_direction =ifelse(Value >0, "Positive", "Negative"),Var =str_to_title(Var)) %>%ggplot(aes(year, Value, group = EPU)) +geom_col(aes(fill = anom_direction), size =0.75, alpha =0.8) +geom_hline(yintercept =0, size =0.5, lty =1) +scale_fill_gmri() +facet_wrap(~Var, ncol =1) +scale_y_continuous(labels = scales::comma_format()) +scale_x_continuous(limits =c(1960, 2020), breaks =seq(1950, 2030, by =10)) +labs(x ="Year", y ="Total Finfish Landings (lb.)",title ="Environmental Drivers") +guides(fill ="none")```#### CPR PCA```{r}#| label: cpr-modes#### CPR Community ##### CPR Plotcpr_indices %>%mutate(anom_direction =ifelse(pca_loading >0, "Positive", "Negative")) %>%ggplot(aes(year, pca_loading, group = pca_mode)) +geom_col(aes(fill = anom_direction), size =0.75, alpha =0.8) +geom_hline(yintercept =0, size =0.5, lty =1) +scale_fill_gmri() +facet_wrap(~pca_mode, ncol =1) +scale_x_continuous(breaks =seq(1950, 2030, by =10)) +labs(x ="Year", y ="PCA Loading", title ="Gulf of Maine CPR PCA") +guides(fill ="none")```:::## Index of WhatOur ability to make statements on community size-structure change in this study is limited by the community that is effectively sampled as part of the groundfish survey, and of species with available weight-at-length information. When looking at the species compositions in each region it is clear that in SNE and MAB a larger proportion of the community biomass was held within a relatively small number of species, primarily the elasmobranchs. Whether or not this is incomplete representation of the community due to the sampling gear, or an accurate representation of the entire community is unclear. In these two regions the biomass spectra information is largely that of the elasmobranch community. Whether or not this is a sufficient fraction of the broader community, what key community members are absent, and what the ecological implications of this are remain unclear.```{r}#| label: proportion-of-fgroup#| fig-height: 8# panel plot the biomass by body sizefgroup_area$weight_bins %>%mutate(survey_area =factor(survey_area, levels = area_levels),strat_lwbio_mill = wtbin_strat_lw_bio /1e6) %>%group_by(Year, survey_area, hare_group) %>%summarise(total_lwbio_mill =sum(strat_lwbio_mill, na.rm = T),.groups ="drop") %>%ggplot(aes(Year, total_lwbio_mill, fill = hare_group)) +geom_col(position ="fill") +facet_wrap(~survey_area, ncol =1) +scale_fill_gmri() +scale_y_continuous(labels =percent_format()) +labs(y ="Fraction of Stratified Biomass", title ="Biomass Allocation by Functional Group",x ="", fill ="Functional Group") +theme(legend.position ="bottom",axis.text.x =element_text(angle =45, hjust =1))```Looking at both the level of noise in the body-size trends and the size spectra parameters, the Gulf of Maine and Georges Bank seem to have the clearest signals coming through. It is likely not coincidental that these two regions also have large proportions of their sampled biomass represented by groundfish and demersal species, functional groups that are best sampled by a bottom trawl survey. Looking specifically at these two regions only, we see two parallel trends of declining body size (length and weight) among the community as a whole. We also simultaneously see a steepening of the size spectra slope. Both of these trends are slightly more severe in the Gulf of Maine.# Supplemental Materials::: panel-tabset## Species Functional Groups```{r}#| label: functional-group-table# Tidy up a key of what is labeled whatspec_key <- catch_size_bins %>%distinct(comname, hare_group) %>%mutate(hare_group =ifelse(is.na(hare_group), "Not Classified", hare_group),hare_group =str_to_title(hare_group),hare_group =fct_drop(hare_group),hare_group =factor(hare_group, levels =c("Coastal", "Demersal", "Elasmobranch", "Groundfish", "Pelagic", "Reef", "Not Classified")) )# Reformat the data to check presence absencecatch_size_bins %>%distinct(comname, survey_area) %>%mutate(Presence ="X") %>%complete(comname, survey_area) %>%left_join(spec_key, by ="comname") %>%mutate(comname =str_to_title(comname),Presence =ifelse(is.na(Presence) ==FALSE, "X", ""),Presence =ifelse(is.na(Presence) ==TRUE, "", Presence)) %>%left_join(area_df, by ="survey_area") %>%select(hare_group, comname, area, Presence) %>%pivot_wider(names_from = area, values_from = Presence) %>%arrange(hare_group, comname) %>%group_by(hare_group) %>%gt() %>% gt::tab_header(title =md("**Functional Group Assignments and Regional Presence/Absence**")) %>%tab_stubhead(label =md("**Functional Group**")) %>%cols_label(comname =md("*Common Name*"),`Gulf of Maine`=md("*Gulf of Maine*"),`Georges Bank`=md("*Georges Bank*"),`Southern New England`=md("*Southern New England*"),`Mid-Atlantic Bight`=md("*Mid-Atlantic Bight*") ) %>%tab_style(style =list(cell_text(weight ="bold")),locations =cells_row_groups() ) %>% gt::tab_source_note(source_note =md("*Functional group assignments adapted from ________*"))```## GARFO Landings Summary```{r}#| label: garfo-landings-summary# Add the labels into the landings data and remove what we don't need there:landings %>%filter(!str_detect(sppname, "CONFIDENTIAL")) %>%mutate(disturbance_era =ifelse(year <=1995, "1964 - 1995", "1996 - 2021"),sppname =str_to_title(sppname)) %>%group_by(survey_area, disturbance_era, sppname) %>%summarise(avg_landings_lb =mean(`landed lbs`, na.rm = T),across(.cols =c(landings_lb =`landed lbs`, value = value), .fns = sum, .names ="total_{.col}")) %>%slice_max(order_by = total_landings_lb, n =5) %>%ungroup() %>%group_by(disturbance_era) %>%gt(rowname_col ="survey_area") %>% gt::tab_header(title =md("**Largest Commercial Fisheries Landings of Northeastern US (by weight)**")) %>%tab_stubhead(label =md("*Harvest Region*")) %>%fmt_number(columns =c(avg_landings_lb, total_landings_lb, total_value),use_seps = T, sep_mark =",",suffixing = T) %>%tab_style(style =list(cell_text(weight ="bold")),locations =cells_row_groups()) %>%cols_label(disturbance_era =md("*Period*"),sppname =md("*Common Name*"),avg_landings_lb =md("*Avg. Annual Landings (lb.)*"),total_landings_lb =md("*Total Landings (lb.)*"),total_value =md("*Total Value ($)*")) %>% gt::tab_source_note(source_note =md("*Landings data obtained from the Greater Atlantic Regional Fishing Office (GARFO)*"))```:::```{r}#| label: ccf-data#| eval: false# Run CCF function to get lagged correlations,# Put that data in one table to plot in a better comparison form# Scale Everything Before CCF,# Add some organization columnsdriver_ccf_prep <- driver_matrix %>%scale() %>%as.data.frame() %>%rownames_to_column(var ="year") %>%pivot_longer(names_to ="spectra_param", values_to ="spectra_values", cols =ends_with("slope") |ends_with("int")) %>%pivot_longer(names_to ="driver_var", values_to ="driver_values", cols =-matches("spectra|year")) %>%mutate(# C. These are the areas associated with the Spectra Featuresxvar_area =case_when(str_detect(spectra_param, "All") ~"All",str_detect(spectra_param, "Georges") ~"Georges Bank",str_detect(spectra_param, "Gulf") ~"Gulf of Maine",str_detect(spectra_param, "Southern") ~"Southern New England",str_detect(spectra_param, "Mid-Atlantic") ~"Mid-Atlantic Bight"),xvar_area =factor(xvar_area, levels =c("All", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")),# D. These are the areas of the driversyvar_area =case_when(str_detect(driver_var, "All") ~"All",str_detect(driver_var, "Georges") ~"Georges Bank",str_detect(driver_var, "Gulf") ~"Gulf of Maine",str_detect(driver_var, "Southern") ~"Southern New England",str_detect(driver_var, "Mid-Atlantic") ~"Mid-Atlantic Bight"),yvar_area =factor(yvar_area, levels =c("All", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight"))) %>%filter(yvar_area == xvar_area |str_detect(yvar_area, "All")) %>%arrange(year, driver_var)# function to grab the correlation data and lag dataget_ccf_vector <-function(x,y){ ccf_data <-ccf(x,y,plot= F)data.frame("acf"= ccf_data$acf,"lag"= ccf_data$lag )}# Walk through each xvariable, and see how it correlates with each yvar# split on xvar# then split on yvars# make sure order is goodccf_relationships <- driver_ccf_prep %>%split(.$spectra_param) %>%map_dfr(function(x_param){ x_param %>%split(.$driver_var) %>%map_dfr(function(driver_y_data){ ccf_df <-get_ccf_vector(x = driver_y_data$spectra_values,y = driver_y_data$driver_values ) }, .id ="driver_var") }, .id ="spectra_param")# Plot themccf_relationships %>%mutate(driver_type =ifelse(str_detect(driver_var, "landings"),"Top-Down","Bottom Up"),param_feature =ifelse(str_detect(spectra_param, "int"),"Intercept","Slope"),region =case_when(str_detect(driver_var, "All") ~"All",str_detect(driver_var, "Georges") ~"Georges Bank",str_detect(driver_var, "Gulf") ~"Gulf of Maine",str_detect(driver_var, "Southern") ~"Southern New England",str_detect(driver_var, "Mid-Atlantic") ~"Mid-Atlantic Bight"),region =factor(region, levels =c("All", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight"))) %>%ggplot(aes(lag, acf, group = driver_var, color = driver_type, fill = driver_type)) +geom_col(alpha =0.5, position ="dodge") +# geom_line(alpha = 0.5, size = 1) +geom_hline(yintercept =0, color ="gray50", size =1) +geom_vline(xintercept =0, color ="black", size =1) +facet_grid(region~param_feature) +scale_x_continuous(limits =c(-12, 12), breaks =seq(from =-12, to =12, by =1)) +theme(legend.position ="bottom") ccf(driver_matrix$Gulf_of_Maine_landings, driver_matrix$Gulf_of_Maine_spectra_slope)```# References